Let $C=C^{\infty}[0,1]$ be the set of functions with derivatives of all orders on $[0,1]$.
Define $n$ linearly independent boundary forms $\{B_j: C\to\mathbb{C}\mid j\in\{1,2,\ldots,n\}\}$ as follows: \begin{align*} B_j\phi &= \sum_{k=0}^{n-1} \left(b_{jk}\phi^{(k)}(0) + \beta_{jk}\phi^{(k)}(1)\right),\quad j\in \{1,2,\ldots, n\}, \end{align*} where $b_{jk}, \beta_{jk}\in\mathbb{R}$ are boundary coefficients. Define \begin{align*} \Phi &= \{\phi\in C: B_j\phi = 0\,\forall j \in \{1,2,\ldots, n\}\} \end{align*} to be the set of $n$ homogeneous boundary conditions \begin{align*} \left(b_{11}\phi^{(0)}(0) + \beta_{11}\phi^{(0)}(1)\right) + \cdots +\left(b_{1n}\phi^{(n-1)}(0) + \beta_{1n}\phi^{(n-1)}(1)\right) &= 0\\ &\vdots\\ \left(b_{n1}\phi^{(0)}(0) + \beta_{n1}\phi^{(0)}(1)\right) + \cdots +\left(b_{nn}\phi^{(n-1)}(0) + \beta_{nn}\phi^{(n-1)}(1)\right) &= 0. \end{align*} Let $S:\Phi\to C$ be the linear differential operator $$S\phi(x) = (-i)^n \phi^{(n)}(x).$$
Let $a\in\mathbb{C}$ be a constant. Consider the homogeneous initial-boundary value problem (IBVP): \begin{alignat*}{2} (\partial_t + aS)q(x,t) &= 0,\quad&(x,t)\in (0,1)\times (0,T)\\ q(x,0)&\in\Phi,\quad &x\in [0,1]\\ q(\cdot, t) &\in\Phi,\quad &t\in [0,T] \end{alignat*} Let $f(x):=q(x,0)$. We reqire that if $n$ is odd then $a=\pm i$ and if $n$ is even then $\Re(a)\geq 0$. We solve the IBVP for $q(x,t)$.
using SymPy
using PyCall
sympy = pyimport("sympy")
# using Roots
using Distributions
# using IntervalArithmetic
# using IntervalRootFinding
using ApproxFun
using Plots
using NBInclude
using QuadGK
import QuadGK.quadgk
# using HCubature
using ApproxFun
using Roots
using Gadfly
using PyPlot
# pygui(true)
using StringLiterals
TOL¶Error tolerance level.
TOL = 1e-05
DIGITS¶The number of digits to display in symbolic expressions.
DIGITS = 3
INFTY¶A number representing infinity in numerical approximations.
INFTY = 10
t¶Free symbol in the unknown function $x(t)$ in the differential equation $Lx=0$.
# t = symbols("t")
sympy_Exprs¶Sample expressions of type addition, multiplication, power, and exponential in SymPy.
x = symbols("x")
sympyAddExpr = 1 + x
sympyMultExpr = 2*x
sympyPowerExpr = x^2
sympyExpExpr = e^x
check_all(array, condition)¶Checks whether all elements in an array satisfy a given condition.
function check_all(array, condition)
for x in array
if !condition(x)
return false
end
end
return true
end
Parameters
array: Arraycondition: Function: BoolReturns
check_all: Booltrue if all elements in array satisfy condition and false if any element does not satisfy the condition.Example
array = [0,1,2,3]
condition = x -> x>2
check_all(array, condition)
check_any(array, condition)¶Checks whether any element in an array satisfy a given condition.
function check_any(array, condition)
for x in array
if condition(x)
return true
end
end
return false
end
Parameters
array: Arraycondition: Function: BoolReturns
check_all: Booltrue if there exists an element in the array that satisfies a given condition and false if no element satisfies the condition.Example
array = [0,1,2,3]
condition = x -> x>2
check_any(array, condition)
set_tol(x, y)¶Sets an appropriate tolerance for checking whether two numbers or two arrays of numbers are approximately equal.
function set_tol(x::Union{Number, Array}, y::Union{Number, Array}; atol = TOL)
if isa(x, Number) && isa(y, Number)
return atol * mean([x y])
elseif isa(x, Array) && isa(y, Array)
if size(x) != size(y)
throw(error("Array dimensions do not match"))
end
# Avoid InexactError() when taking norm()
x = convert(Array{Complex}, x)
y = convert(Array{Complex}, y)
return atol * (norm(x,2) + norm(y,2))
else
throw(error("Invalid input"))
end
end
Parameters
x, y: Number or Array of NumberReturns
set_tol: Numberx and y are considered approximately equal (entry-wise if x, y are arrays).Example
x = 14
y = 21
println("set_tol($x, $y) = $(set_tol(x,y))")
A = [4 0.6; 3 2]
B = [5 1; 10 3]
set_tol(A, B)
println("set_tol($A, $B) = $(set_tol(A,B))")
evaluate(func, a)¶Evaluate a univariate function or an array of them at a given value.
function evaluate(func::Union{Function, SymPy.Sym, Number}, a::Number)
if isa(func, Function)
funcA = func(a)
elseif isa(func, SymPy.Sym) # SymPy.Sym must come before Number because SymPy.Sym will be recognized as Number
freeSymbols = free_symbols(func)
if length(freeSymbols) > 1
throw(error("func should be univariate"))
elseif length(freeSymbols) == 1
t = free_symbols(func)[1,1]
if isa(a, SymPy.Sym) # if x is SymPy.Sym, do not convert result to Number to preserve pretty printing
funcA = subs(func, t, a)
else
funcA = SymPy.N(subs(func, t, a))
end
else
funcA = func
end
else # func is Number
funcA = func
end
return funcA
end
Parameters
func: Function, SymPy.Sym, Number, or ArraySymPy.Sym is absorbed into Number.x: Numberfunc is to be evaluated at.t: SymPy.Sym func if func is a SymPy.Sym object or an array of them.Returns
evaluate: Numberfunc at x.Example
func = x -> x+1
xVal = 2
println("func($xVal) = $(evaluate(func, xVal))")
x = symbols("x")
func = x+1
println("func($xVal) = $(evaluate(func, xVal))")
a = symbols("a")
println("func($a) = $(evaluate(func, a))")
x = symbols("x")
array = [2x 1; x^3 x]
a = symbols("a")
evaluate.(array, a)
partition(n)¶Generate ordered two-integer partitions ($0$ included) of a given non-negative integer.
function partition(n::Int)
if n < 0
throw(error("Non-negative n required"))
end
output = []
for i = 0:n
j = n - i
push!(output, (i,j))
end
return output
end
Parameters
n: IntReturns
partition: Array of Tuple{Int, Int}n. Note that a tuple is ordered, and (i, j) and (j, i) are considered distinct partitions.Example
n = 5
partition(5)
get_deriv(u, k)¶Constructs the symbolic expression for the $k$th derivative of a univariate function u.
function get_deriv(u::Union{SymPy.Sym, Number}, k::Int)
if k < 0
throw(error("Non-negative k required"))
end
if isa(u, SymPy.Sym)
freeSymbols = free_symbols(u)
if length(freeSymbols) > 1
throw(error("u should be univariate"))
elseif length(freeSymbols) == 1
t = freeSymbols[1]
y = u
for i = 1:k
newY = diff(y, t)
y = newY
end
return y
else
if k > 0
return 0
else
return u
end
end
else
if k > 0
return 0
else
return u
end
end
end
Parameters
u: SymPy.Sym or Numberk: IntReturns
get_deriv: SymPy.Symu.Example
u = 3
k = 1
println("get_deriv($u, $k) = $(get_deriv(u, k))")
t = symbols("t")
u = t^2+2t+3
get_deriv(u, k)
add_func(f, g)¶Computes the sum of two functions using the function addition given by $(f + g)(x) := f(x) + g(x)$.
function add_func(f::Union{Number, Function}, g::Union{Number, Function})
function h(x)
if isa(f, Number)
if isa(g, Number)
return f + g
else
return f + g(x)
end
elseif isa(f, Function)
if isa(g, Number)
return f(x) + g
else
return f(x) + g(x)
end
end
end
return h
end
Parameters
f, g: Function or NumberReturns
add_func: Function or Numberf and g.Example
f(x) = x^3+1
g(x) = 4x
x = 2
add_func(f, g)(x) == f(x) + g(x)
mult_func(f, g)¶Computes the product of two functions using the function multiplication given by $(f \cdot g)(x) := f(x) \cdot g(x)$.
function mult_func(f::Union{Number, Function}, g::Union{Number, Function})
function h(x)
if isa(f, Number)
if isa(g, Number)
return f * g
else
return f * g(x)
end
elseif isa(f, Function)
if isa(g, Number)
return f(x) * g
else
return f(x) * g(x)
end
end
end
return h
end
Parameters
f, g: Function or NumberReturns
add_func: Function or Numberf and g.Example
f(x) = x^3+1
g(x) = 4x
x = 2
mult_func(f, g)(x) == f(x) * g(x)
sym_to_func(sym)¶Converts a univariate symbolic expression to a function.
function sym_to_func(sym::Union{SymPy.Sym, Number})
try
freeSymbols = free_symbols(sym)
if length(freeSymbols) > 1
throw(error("sym should be univariate"))
else
function func(x)
if length(freeSymbols) == 0
result = SymPy.N(sym)
else
result = SymPy.N(subs(sym, freeSymbols[1], x))
end
return result
end
return func
end
catch
function func(x)
return sym
end
return func
end
end
Parameters
sym: SymPy.Sym or NumberReturns
sym_to_func: Functionsym.Example
t = symbols("t")
sym = t^2+t
tVal = 2
println("sym_to_func($sym)($tVal) = $(sym_to_func(sym)(tVal))")
sym = [1 t; t^2 t^3]
sym_to_func.(sym)[2,1](2)
println("sym_to_func($sym[2,1])($tVal) = $(sym_to_func.(sym[2,1])(tVal))")
sym = 2
println("sym_to_func($sym)($tVal) = $(sym_to_func(sym)(tVal))")
prettyRound(x, digits)¶Round a number at a given number of digits.
function prettyRound(x::Number; digits::Int = DIGITS)
if isa(x, Int)
return x
elseif isa(x, Real)
if isa(x, Rational) || isa(x, Irrational) # If x is rational or irrational numbers like e, pi
return x
elseif round(abs(x), digits) == floor(abs(x))
return Int(round(x))
else
return round(x, digits)
# return rationalize(x)
end
elseif isa(x, Complex)
roundedReal = prettyRound(real(x), digits = digits)
roundedComplex = prettyRound(imag(x), digits = digits)
return roundedReal + im*roundedComplex
else
return round(x, digits)
end
end
Parameters
x: Numberdigits*: IntDIGITS.Returns
prettyRound: Numberx rounded to the digits decimal places.Example
# x = 0.0000001*im
x = 4.854572304702339 - 49.023298458074294im
prettyRound(x)
println("prettyRound($x) = $(prettyRound(x))")
x = [1/3 1/4 0]
println("prettyRound($x) = $(prettyRound.(x))")
x = [1//3 1//4 0//1]
println("prettyRound($x) = $(prettyRound.(x))")
prettyPrint(x)¶Prints a symbolic scalar pretty-rounded floating numbers.
function prettyPrint(x::Union{Number, SymPy.Sym})
expr = x
if isa(expr, SymPy.Sym)
prettyExpr = expr
for a in sympy[:preorder_traversal](expr)
if length(free_symbols(a)) == 0 && length(args(a)) == 0
if !(a in [e, PI]) && length(intersect(args(a), [e, PI])) == 0 # keep the transcendental numbers as symbols
prettyA = prettyRound.(SymPy.N(a))
prettyExpr = subs(prettyExpr, (a, prettyA))
end
end
end
else
prettyExpr = prettyRound.(expr)
prettyExpr = convert(SymPy.Sym, prettyExpr)
end
return prettyExpr
end
Parameters
x: Number or SymPy.SymReturns
prettyPrint: SymPy.Symx with pretty-rounded floating numbers.Example
t = symbols("t")
x = 1//3*t + e^(2t) + PI/2 + 1.00001*im + sympy[:sqrt](2) + 1//6
prettyPrint(x)
t = symbols("t")
x = [1//3*t PI; 1//6*t^2+t 1]
prettyPrint.(x)
x = [0.0+1.0*im 1.0+0.0*im -1.0+0.0*im 0.0-1.0*im]
prettyPrint.(x)
is_approxLess(x, y; atol)¶Checks whether $x$ is approximately less than $y$ within a tolerance. That is, whether $x\not\approx y$ and $x<y$.
function is_approxLess(x::Number, y::Number; atol = TOL)
return !isapprox(x,y; atol = atol) && x < y
end
Parameters
x, y: Numberatol*: NumberTOL.Returns
is_approxLess: Booltrue if $x\not\approx y$ within atol and $x<y$, and false otherwise.Example
x = 1
y = x + TOL/2
println("is_approxLess($x,$y) = $(is_approxLess(x,y))")
y = x + TOL*2
println("is_approxLess($x,$y) = $(is_approxLess(x,y))")
is_approx(x, y; atol)¶Checks whether $x$ is approximately equal to $y$ within a tolerance.
function is_approx(x::Number, y::Number; atol = TOL)
return isapprox(x, y; atol = atol)
end
Parameters
x, y: Numberatol*: NumberTOL.Returns
is_approx: Booltrue if $x\approx y$ within atol and false otherwise.Example
x = 1
y = x + TOL/2
println("is_approx($x,$y) = $(is_approx(x,y))")
y = x + TOL*2
println("is_approx($x,$y) = $(is_approx(x,y))")
argument(z)¶Finds the argument of a complex number in $[0,2\pi)$.
function argument(z::Number)
if angle(z) >= 0 # in [0,pi]
return angle(z)
else
# angle(z) is in (-pi, 0]
# Shift it up to (pi,2pi]
argument = 2pi + angle(z) # This is in (pi,2pi]
if is_approx(argument, 2pi) # Map 2pi to 0
return 0
else
return argument # This is now in [0,2pi)
end
end
end
Parameters
z: NumberReturns
argument: NumberExample
z = 1
println("argument($z) = $(argument(z))")
z = -1-im
println("argument($z) = $(argument(z))")
trace_contour(a, n, infty, sampleSize)¶Plots the contour sectors encircled by $\Gamma_a^+$, $\Gamma_a^-$, $\Gamma_0^+$, and $\Gamma_0^-$ defined as \begin{align*} \Gamma_a^{\pm} &:= \partial(\{\lambda\in\mathbb{C}^{\pm}:\, \Re(a\lambda^n)>0\}\setminus \bigcup_{\substack{\sigma\in\mathbb{C};\\\Delta(\sigma)=0}} D(\sigma, 2\epsilon))\quad\mbox{($D$ for disk)},\\ \Gamma_0^+ &:= \bigcup_{\substack{\sigma\in\overline{\mathbb{C}^+};\\ \Delta(\sigma)=0}} C(\sigma, \epsilon)\quad\mbox{($C$ for circle)},\\ \Gamma_0^- &:= \bigcup_{\substack{\sigma\in\mathbb{C}^-;\\ \Delta(\sigma)=0}} C(\sigma, \epsilon) \end{align*} by sampling points.
function trace_contour(a::Number, n::Int, sampleSize::Int; infty = INFTY)
lambdaVec = []
for counter = 1:sampleSize
x = rand(Uniform(-infty,infty), 1, 1)[1]
y = rand(Uniform(-infty,infty), 1, 1)[1]
lambda = x + y*im
if real(a*lambda^n)>0
append!(lambdaVec, lambda)
end
end
Gadfly.plot(x=real(lambdaVec), y=imag(lambdaVec), Guide.xlabel("Re"), Guide.ylabel("Im"), Coord.Cartesian(ymin=-infty,ymax=infty, xmin=-infty, xmax=infty, fixed = true))
end
Parameters
a: Numbern: Intinfty: NumberINFTY.sampleSize: IntReturns
trace_contour: NoneExample
n = 2
a = 1
sampleSize = 1000
trace_contour(a, n, sampleSize)
n = 3
a = im
sampleSize = 1000
trace_contour(a, n, sampleSize)
n = 3
a = -im
sampleSize = 1000
trace_contour(a, n, sampleSize)
get_distancePointLine(z, theta)¶Finds the distance between a complex number and a line through the origin given by an angle in $[0,2pi)$.
function get_distancePointLine(z::Number, theta::Number)
if theta >= 2pi && theta < 0
throw(error("Theta must be in [0,2pi)"))
else
if is_approx(argument(z), theta)
return 0
else
x0, y0 = real(z), imag(z)
if is_approx(theta, pi/2) || is_approx(theta, 3pi/2)
return abs(x0)
elseif is_approx(theta, 0) || is_approx(theta, 2pi)
return abs(y0)
else
k = tan(theta)
x = (y0+1/k*x0)/(k+1/k)
y = k*x
distance = norm(z-(x+im*y))
return distance
end
end
end
end
Parameters
z: Numbertheta: NumberReturns
get_distancePointLine: Numberz and the line characterized by theta.Example
z = 1
theta = pi/3
println("get_distancePointLine($z, $theta) = $(get_distancePointLine(z, theta))") # sqrt(3)/2
StructDefinitionError¶A struct definition error type is the class of all errors in struct definitions.
struct StructDefinitionError <: Exception
msg::String
end
SymLinearDifferentialOperator(symPFunctions, interval, t)¶A symbolic linear differential operator of order $n$ is encoded by an $1 \times (n+1)$ array of symbolic expressions with at most one free symbol, an interval $[a,b]$, and that free symbol.
struct SymLinearDifferentialOperator
# Entries in the array should be SymPy.Sym or Number. SymPy.Sym seems to be a subtype of Number, i.e., Array{Union{Number,SymPy.Sym}} returns Array{Number}. But specifying symPFunctions as Array{Number,2} gives a MethodError when the entries are Sympy.Sym objects.
symPFunctions::Array
interval::Tuple{Number,Number}
t::SymPy.Sym
SymLinearDifferentialOperator(symPFunctions::Array, interval::Tuple{Number,Number}, t::SymPy.Sym) =
try
symL = new(symPFunctions, interval, t)
check_symLinearDifferentialOperator_input(symL)
return symL
catch err
throw(err)
end
end
function check_symLinearDifferentialOperator_input(symL::SymLinearDifferentialOperator)
symPFunctions, (a,b), t = symL.symPFunctions, symL.interval, symL.t
for symPFunc in symPFunctions
if isa(symPFunc, SymPy.Sym)
if size(free_symbols(symPFunc)) != (1,) && size(free_symbols(symPFunc)) != (0,)
throw(StructDefinitionError(:"Only one free symbol is allowed in symP_k"))
end
elseif !isa(symPFunc, Number)
throw(StructDefinitionError(:"symP_k should be SymPy.Sym or Number"))
end
end
return true
end
Parameters
symPFunctions: Array of SymPy.Sym or Numberinterval: Tuple{Number, Number} t: SymPy.Sym symPFunctions.Returns
SymLinearDifferentialOperatorSymLinearDifferentialOperator of order $n$ with attributes symPFunctions, interval, and t.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
LinearDifferentialOperator(pFunctions, interval, symL)¶A linear differential operator $L$ of order $n$ given by $$Lx = p_0x^{(n)} + p_1x^{(n-1)} + \cdots + p_{n-1}x^{(1)} + p_n x$$ is encoded by an $1 \times (n+1)$ array of univariate functions, an interval $[a,b]$, and its symbolic expression.
# symL is an attribute of L that needs to be input by the user. There are checks to make sure symL is indeed the symbolic version of L.
# Principle: Functionalities of Julia Functions >= Functionalities of SymPy. If p_k has no SymPy representation, the only consequence should be that outputs by functions that take L as arugment has no symbolic expression. E.g., we allow L.pFunctions and L.symL.pFunctions to differ.
struct LinearDifferentialOperator
pFunctions::Array # Array of julia functions or numbers representing constant functions
interval::Tuple{Number,Number}
symL::SymLinearDifferentialOperator
LinearDifferentialOperator(pFunctions::Array, interval::Tuple{Number,Number}, symL::SymLinearDifferentialOperator) =
try
L = new(pFunctions, interval, symL)
check_linearDifferentialOperator_input(L)
return L
catch err
throw(err)
end
end
# Assume symFunc has only one free symbol, as required by the definition of SymLinearDifferentialOperator.
# That is, assume the input symFunc comes from SymLinearDifferentialOperator.
function check_func_sym_equal(func::Union{Function,Number}, symFunc, interval::Tuple{Number,Number}, t::SymPy.Sym) # symFunc should be Union{SymPy.Sym, Number}, but somehow SymPy.Sym gets ignored
(a,b) = interval
# Randomly sample 1000 points from (a,b) and check if func and symFunc agree on them
for i = 1:1000
# Check endpoints
if i == 1
x = a
elseif i == 2
x = b
else
x = rand(Uniform(a,b), 1)[1,1]
end
funcEvalX = evaluate.(func, x)
if isa(symFunc, SymPy.Sym)
symFuncEvalX = SymPy.N(subs(symFunc,t,x))
# N() converts SymPy.Sym to Number
# https://docs.sympy.org/latest/modules/evalf.html
# subs() works no matter symFunc is Number or SymPy.Sym
else
symFuncEvalX = symFunc
end
tol = set_tol(funcEvalX, symFuncEvalX)
if !isapprox(real(funcEvalX), real(symFuncEvalX); atol = real(tol)) ||
!isapprox(imag(funcEvalX), imag(symFuncEvalX); atol = imag(tol))
println("x = $x")
println("symFunc = $symFunc")
println("funcEvalX = $funcEvalX")
println("symFuncEvalX = $symFuncEvalX")
return false
end
end
return true
end
# Check whether the inputs of L are valid.
function check_linearDifferentialOperator_input(L::LinearDifferentialOperator)
pFunctions, (a,b), symL = L.pFunctions, L.interval, L.symL
symPFunctions, t = symL.symPFunctions, symL.t
# domainC = Complex(a..b, 0..0) # Domain [a,b] represented in the complex plane
p0 = pFunctions[1]
# p0Chebyshev = Fun(p0, a..b) # Chebysev polynomial approximation of p0 on [a,b]
if !check_all(pFunctions, pFunc -> (isa(pFunc, Function) || isa(pFunc, Number)))
throw(StructDefinitionError(:"p_k should be Function or Number"))
elseif length(pFunctions) != length(symPFunctions)
throw(StructDefinitionError(:"Number of p_k and symP_k do not match"))
elseif (a,b) != symL.interval
throw(StructDefinitionError(:"Intervals of L and symL do not match"))
# # Assume p_k are in C^{n-k}. Check whether p0 vanishes on [a,b].
# # roots() in IntervalRootFinding doesn't work if p0 is sth like t*im - 2*im. Neither does find_zero() in Roots.
# # ApproxFun.roots()
# elseif (isa(p0, Function) && (!isempty(roots(p0Chebyshev)) || all(x->x>b, roots(p0Chebyshev)) || all(x->x<b, roots(p0Chebyshev)) || p0(a) == 0 || p0(b) == 0)) || p0 == 0
# throw(StructDefinitionError(:"p0 vanishes on [a,b]"))
elseif !all(i -> check_func_sym_equal(pFunctions[i], symPFunctions[i], (a,b), t), 1:length(pFunctions))
# throw(StructDefinitionError(:"symP_k does not agree with p_k on [a,b]"))
warn("symP_k does not agree with p_k on [a,b]") # Make this a warning instead of an error because the functionalities of Julia Functions may be more than those of SymPy objects; we do not want to compromise the functionalities of LinearDifferentialOperator because of the restrictions on SymPy.
else
return true
end
end
Parameters
pFunctions: Array of Function or Number interval: Tuple{Number, Number}(a, b) corresponding to the real interval $[a,b]$ on which the differential operator $L$ is defined.symL: SymLinearDifferentialOperatorReturns
LinearDifferentialOperatorLinearDifferentialOperator of order $n$ with attributes pFunctions, interval, and symL.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
# Direct construction
pFunctions = [t->1 t->t+1 t->t^2+t+1]
L = LinearDifferentialOperator(pFunctions, interval, symL)
VectorBoundaryForm(M, N)¶A set of homogeneous boundary conditions in vector form $$Ux = \begin{bmatrix}U_1\\\vdots\\ U_m\end{bmatrix}x = \begin{bmatrix}\sum_{j=1}^n M_{1j}x^{(j-1)}(a) + N_{1j}x^{(j-1)}(b)\\\vdots\\ \sum_{j=1}^n M_{mj}x^{(j-1)}(a) + N_{mj}x^{(j-1)}(b)\end{bmatrix} = \begin{bmatrix}0\\\vdots\\ 0\end{bmatrix}$$ is encoded by an ordered pair of two linearly independent $m\times n$ matrices $(M, N)$ where $$M = \begin{bmatrix}M_{11} & \cdots & M_{1n}\\ \vdots & \ddots & \vdots\\ M_{m1} & \cdots & M_{mn}\end{bmatrix},\quad N = \begin{bmatrix}N_{11} & \cdots & N_{1n}\\ \vdots & \ddots & \vdots\\ N_{m1} & \cdots & N_{mn}\end{bmatrix}.$$
struct VectorBoundaryForm
M::Array # Why can't I specify Array{Number,2} without having a MethodError?
N::Array
VectorBoundaryForm(M::Array, N::Array) =
try
U = new(M, N)
check_vectorBoundaryForm_input(U)
return U
catch err
throw(err)
end
end
# Check whether the input matrices that characterize U are valid
function check_vectorBoundaryForm_input(U::VectorBoundaryForm)
# M, N = U.M, U.N
# Avoid Inexact() error when taking rank()
M = convert(Array{Complex}, U.M)
N = convert(Array{Complex}, U.N)
if !(check_all(U.M, x -> isa(x, Number)) && check_all(U.N, x -> isa(x, Number)))
throw(StructDefinitionError(:"Entries of M, N should be Number"))
elseif size(U.M) != size(U.N)
throw(StructDefinitionError(:"M, N dimensions do not match"))
elseif size(U.M)[1] != size(U.M)[2]
throw(StructDefinitionError(:"M, N should be square matrices"))
elseif rank(hcat(M, N)) != size(M)[1] # rank() throws weird "InexactError()" when taking some complex matrices
throw(StructDefinitionError(:"Boundary operators not linearly independent"))
else
return true
end
end
Parameters
M, N: Array of NumberReturns
VectorBoundaryFormVectorBoundaryForm with attributes M and N.Example
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
Constructs a valid adjoint boundary condition from a given (homogeneous) boundary condition based on Chapter 11 in Theory of Ordinary Differential Equations (Coddington & Levinson).
get_L(symL)¶Constructs a LinearDifferentialOperator from a given SymLinearDifferentialOperator.
function get_L(symL::SymLinearDifferentialOperator)
symPFunctions, (a,b), t = symL.symPFunctions, symL.interval, symL.t
if check_all(symPFunctions, x->!isa(x, SymPy.Sym))
pFunctions = symPFunctions
else
pFunctions = sym_to_func.(symPFunctions)
end
L = LinearDifferentialOperator(pFunctions, (a,b), symL)
return L
end
Parameters
symL: SymLinearDifferentialOperatorReturns
get_L: LinearDifferentialOperatorsymL.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
get_URank(U)¶Computes the rank of a vector boundary form $U$ by computing the equivalent $\text{rank}(M:N)$, where $M, N$ are the matrices associated with $U$ and $$(M:N) = \begin{bmatrix}M_{11} & \cdots & M_{1n} & N_{11} & \cdots & N_{1n}\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots\\ M_{m1} & \cdots & M_{mn} & N_{m1} & \cdots & N_{mn}\end{bmatrix}.$$
function get_URank(U::VectorBoundaryForm)
# Avoid InexactError() when taking hcat() and rank()
M = convert(Array{Complex}, U.M)
N = convert(Array{Complex}, U.N)
MHcatN = hcat(M, N)
return rank(MHcatN)
end
Parameters
U: VectoBoundaryFormReturns
get_URank: NumberU.Example
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
get_URank(U)
get_Uc(U)¶Given vector boundary form $U = \begin{bmatrix}U_1\\ \vdots\\ U_m\end{bmatrix}$ of rank $m$, finds a complementary form $U_c = \begin{bmatrix}U_{m+1}\\ \vdots\\ U_{2n}\end{bmatrix}$ of rank $2n-m$ such that $\begin{bmatrix}U_1\\ \vdots\\ U_{2n}\end{bmatrix}$ has rank $2n$.
function get_Uc(U::VectorBoundaryForm)
try
check_vectorBoundaryForm_input(U)
n = get_URank(U)
I = complex(eye(2*n))
M, N = U.M, U.N
MHcatN = hcat(M, N)
# Avoid InexactError() when taking rank()
mat = convert(Array{Complex}, MHcatN)
for i = 1:(2*n)
newMat = vcat(mat, I[i:i,:])
newMat = convert(Array{Complex}, newMat)
if rank(newMat) == rank(mat) + 1
mat = newMat
end
end
UcHcat = mat[(n+1):(2n),:]
Uc = VectorBoundaryForm(UcHcat[:,1:n], UcHcat[:,(n+1):(2n)])
return Uc
catch err
return err
end
end
Parameters
U: VectoBoundaryFormReturns
get_Uc: VectorBoundaryFormU.Example
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
Uc = get_Uc(U)
get_H(U, Uc)¶Given a vector boundary form $U$ and a complementary vector boundary form $U_c$, constructs $$H = \begin{bmatrix}M&N\\ M_c & N_c\end{bmatrix},$$ where $M, N$ are the matrices associated with $U$ and $M_c, N_c$ are associated with $U_c$.
function get_H(U::VectorBoundaryForm, Uc::VectorBoundaryForm)
MHcatN = hcat(convert(Array{Complex}, U.M), convert(Array{Complex}, U.N))
McHcatNc = hcat(convert(Array{Complex}, Uc.M), convert(Array{Complex}, Uc.N))
H = vcat(MHcatN, McHcatNc)
return H
end
Parameters
U: VectorBoundaryFormUc: VectorBoundaryFormU.Returns
get_H: ArrayExample
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
Uc = get_Uc(U)
get_H(U, Uc)
get_pDerivMatrix(L; symbolic)¶Given a LinearDifferentialOperator L where L.pFunctions is the array
$$[p_0, p_1, \ldots, p_n],$$
constructs an $n\times n$ matrix whose $(i+1)(j+1)$-entry is a function corresponding to the $j$th derivative of $p_i$:
$$\begin{bmatrix}p_0 & \cdots & p_0^{(n-1)}\\ \vdots & \ddots & \vdots\\ p_{n-1} & \cdots & p_{n-1}^{(n-1)}\end{bmatrix}.$$
function get_pDerivMatrix(L::LinearDifferentialOperator; symbolic = false)
if symbolic
symL = L.symL
symPFunctions, t = symL.symPFunctions, symL.t
n = length(symPFunctions)-1
symPDerivMatrix = Array{SymPy.Sym}(n,n)
pFunctionSymbols = symPFunctions
for i in 0:(n-1)
for j in 0:(n-1)
index, degree = i, j
symP = pFunctionSymbols[index+1]
# If symP is not a Sympy.Sym object (e.g., is a Number instead), then cannot use get_deriv()
if !isa(symP, SymPy.Sym)
if degree > 0
symPDeriv = 0
else
symPDeriv = symP
end
else
symPDeriv = get_deriv(symP, degree)
end
symPDerivMatrix[i+1,j+1] = symPDeriv
end
end
return symPDerivMatrix
else
symPDerivMatrix = get_pDerivMatrix(L; symbolic = true)
n = length(L.pFunctions)-1
pDerivMatrix = sym_to_func.(symPDerivMatrix)
end
return pDerivMatrix
end
Parameters
L: LinearDifferentialOperatorpDerivMatrix is to be constructed.symbolic*: BoolReturns
get_pDerivMatrix: Array of Function, Number or SymPy.SymL.pFunctions[i]) if symbolic = false, orL.symL.symPFunctions[i]) if symbolic = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
# pFunctions = [t->1 t->t+1 t->t^2+t+1]
# L = LinearDifferentialOperator(pFunctions, interval, symL)
L = get_L(symL)
tVal = 2
pDerivMatrix = get_pDerivMatrix(L; symbolic = false)
println("pDerivMatrix($tVal) = $(evaluate.(pDerivMatrix,tVal))")
symPDerivMatrix = get_pDerivMatrix(L; symbolic = true)
get_Bjk(L, j, k; symbolic, pDerivMatrix)¶Given a LinearDifferentialOperator L of order $n$, for $j, k \in \{1,\ldots,n\}$, computes $B_{jk}$ defined as
$$B_{jk}(t) := \sum_{\ell=j-1}^{n-k}\binom{\ell}{j-1}p^{(\ell-j+1)}_{n-k-\ell}(t)(-1)^\ell.$$
function get_Bjk(L::LinearDifferentialOperator, j::Int, k::Int; symbolic = false, pDerivMatrix = get_pDerivMatrix(L; symbolic = symbolic))
n = length(L.pFunctions)-1
if j <= 0 || j > n || k <= 0 || k > n
throw("j, k should be in {1, ..., n}")
end
sum = 0
if symbolic
symPDerivMatrix = get_pDerivMatrix(L; symbolic = true)
for l = (j-1):(n-k)
summand = binomial(l, j-1) * symPDerivMatrix[n-k-l+1, l-j+1+1] * (-1)^l
sum += summand
end
else
for l = (j-1):(n-k)
summand = mult_func(binomial(l, j-1) * (-1)^l, pDerivMatrix[n-k-l+1, l-j+1+1])
sum = add_func(sum, summand)
end
end
return sum
end
Parameters
L: LinearDifferentialOperatorL.pFunctions are to become the $p_{n-k-l}^{l-j+1}$ in $B_{jk}(t)$.j, k: Intsymbolic*: BoolpDerivMatrix*: Arraysymbolic = false, an $n\times n$ matrix whose $(i+1)(j+1)$-entry is the $j$th derivative of $p_i$ (L.pFunctions[i]) implemented as a Function, Number, or SymPy.Sym. Default to the output of get_pDerivMatrix(L).Returns
get_Bjk: SymPy.Sym, Function, or NumberFunction if symbolic = false, orSymPy.Sym object if symbolic = true, where each $p_i$ is the generic expression $p_i(t)$.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
# pFunctions = [t->1 t->t+1 t->t^2+t+1]
# L = LinearDifferentialOperator(pFunctions, interval, symL)
L = get_L(symL)
# pDerivMatrix = [1 0; t->1+t t->1]
j, k = 1, 1
BjkSym = get_Bjk(L, j, k; symbolic = true)
tVal = 1
Bjk = get_Bjk(L, j, k; symbolic = false)
println("Bjk($tVal) = $(Bjk(tVal))")
println("BjkSym($tVal) = $(evaluate.(BjkSym, tVal))")
get_B(L; symbolic, pDerivMatrix)¶Given a LinearDifferentialOperator L where L.pFunctions is the array
$$[p_0, p_1, \ldots, p_n],$$
constructs the matrix $B(t)$ whose $ij$-entry is given by
$$B_{jk}(t) := \sum_{\ell=j-1}^{n-k}\binom{\ell}{j-1}p^{(\ell-j+1)}_{n-k-\ell}(t)(-1)^\ell.$$
function get_B(L::LinearDifferentialOperator; symbolic = false, pDerivMatrix = get_pDerivMatrix(L; symbolic = symbolic))
n = length(L.pFunctions)-1
B = Array{Union{Function, Number, SymPy.Sym}}(n,n)
for j = 1:n
for k = 1:n
B[j,k] = get_Bjk(L, j, k; symbolic = symbolic, pDerivMatrix = pDerivMatrix)
end
end
return B
end
Parameters
L: LinearDifferentialOperatorL.pFunctions are to become the $p_{n-k-l}^{l-j+1}$ in $B_{jk}(t)$.symbolic*: Boolsubstitute*: Boolsymbolic = true, boolean indicating whether to substitute the symbolic expression of $p_i$ in L.pFunctions for the generic expression $p_i(t)$ created using SymFunction("pi")(t). If symbolic = false, the value of substitute does not matter.pDerivMatrix*: Arraysymbolic = false, the non-symbolic version of symPDerivMatrix, i.e., an $n\times n$ matrix whose $(i+1)(j+1)$-entry is the $j$th derivative of $p_i$ (L.pFunctions[i]) implemented as a Function or Number.Returns
get_B: Array of Function, SymPy.Sym, or NumberFunction if symbolic = false, orSymPy.Sym object if symbolic = true, where each $p_i$ issubstitute = false, orL.symL.symPFunctions[i]) if substitute = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
BSym = get_B(L; symbolic = true)
# Only Array{SymPy.Sym} would be automatically pretty-printed
# Enfore pretty-printing
prettyPrint.(BSym)
B = get_B(L; symbolic = false)
tVal = 1
println("B($tVal) = $(evaluate.(B, tVal))")
println("BSym($tVal) = $(evaluate.(BSym, tVal))")
get_BHat(L, B)¶Given a LinearDifferentialOperator L where L.pFunctions is the array
$$[p_0, p_1, \ldots, p_n]$$
and L.interval is $[a,b]$, constructs $\hat{B}$ defined as the block matrix
$$\hat{B}:=\begin{bmatrix}-B(a) & 0_n\\0_n & B(b)\end{bmatrix}.$$
function get_BHat(L::LinearDifferentialOperator, B::Array)
# if check_any(B, x->isa(x, SymPy.Sym))
# throw("Entries of B should be Function or Number")
# end
pFunctions, (a,b) = L.pFunctions, L.interval
n = length(pFunctions)-1
BHat = Array{Complex}(2n,2n)
BEvalA = evaluate.(B, a)
BEvalB = evaluate.(B, b)
BHat[1:n,1:n] = -BEvalA
BHat[(n+1):(2n),(n+1):(2n)] = BEvalB
BHat[1:n, (n+1):(2n)] = 0
BHat[(n+1):(2n), 1:n] = 0
return BHat
end
Parameters
L: LinearDifferentialOperatorL.pFunctions are to become the $p_{n-k-l}^{l-j+1}$ in $B_{jk}(t)$.B: Array of Numberget_B(L; symbolic = false).Returns
get_BHat: Array of NumberExample
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
B = get_B(L; symbolic = false)
get_BHat(L, B)
get_J(BHat, H)¶Given $\hat{B}$ and $H$, constructs $J$ defined as $$J:=(\hat{B}H^{-1})^\star$$ where $^*$ denotes conjugate transpose.
function get_J(BHat, H)
n = size(H)[1]
H = convert(Array{Complex}, H)
J = (BHat * inv(H))'
# J = convert(Array{Complex}, J)
return J
end
Parameters
BHat: Arrayget_BHat().H: Arrayget_H().Returns
get_J: ArrayExample
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
B = get_B(L; symbolic = false)
BHat = get_BHat(L, B)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
Uc = get_Uc(U)
H = get_H(U, Uc)
get_J(BHat, H)
get_adjointCand(J)¶Given $J$, constructs a candidate adjoint vector boundary form $U^+$ from two matrices $P^\star$, $Q^\star$, which are the lower-left $n\times n$ submatrix of $J$, and the lower-right $n\times n$ submatrix of $J$, respectively.
function get_adjointCand(J)
n = convert(Int, size(J)[1]/2)
J = convert(Array{Complex}, J)
PStar = J[(n+1):2n,1:n]
QStar = J[(n+1):2n, (n+1):2n]
adjointU = VectorBoundaryForm(PStar, QStar)
return adjointU
end
Parameters
J: Arrayget_J.Returns
get_adjoint: VectorBoundaryFormExample
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
B = get_B(L; symbolic = false)
BHat = get_BHat(L, B)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
Uc = get_Uc(U)
H = get_H(U, Uc)
J = get_J(BHat, H)
adjoint = get_adjointCand(J)
get_xi(L; symbolic, xSym)¶Given a LinearDifferentialOperator L of order $n$ in the differential equation $Lx=0$, constructs $\xi(t)$, which is defined as the vector of derivatives of $x(t)$
$$\xi(t) := \begin{bmatrix}x(t)\\ x^{(1)}(t)\\ x^{(2)}(t)\\ \vdots\\ x^{(n-1)}(t)\end{bmatrix}.$$
function get_xi(L::LinearDifferentialOperator; symbolic = true, xSym= nothing)
if symbolic
t = L.symL.t
n = length(L.pFunctions)-1
symXi = Array{SymPy.Sym}(n,1)
if isa(xSym, Void)
throw(error("xSymrequired"))
else
for i = 1:n
symXi[i] = get_deriv(xSym, i-1)
end
return symXi
end
else
if isa(xSym, Void)
throw(error("xSym required"))
elseif !isa(xSym, SymPy.Sym) && !isa(xSym, Number)
throw(error("xSym should be SymPy.Sym or Number"))
else
symXi = get_xi(L; symbolic = true, xSym = xSym)
xi = sym_to_func.(symXi)
end
end
end
Parameters
L: LinearDifferentialOperatorsymbolic: Boolsubstitute*: Boolsymbolic = true, boolean indicating whether to substitute the symbolic expression of $x(t)$ for the generic expression created using SymFunction.xSym*: SymPy.Symsubstitute = true, symbolic expression of $x(t)$ to replace the generic expression with.Returns
get_xi: Array of SymPy.Symsubstitute = false, orsubstitute = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
xSym = t^2+2t
symXi = get_xi(L; symbolic = true, xSym = xSym)
xi = get_xi(L; symbolic=false, xSym = xSym)
tVal = 1
println("xi($tVal) = $(evaluate.(xi, tVal))")
println("symXi($tVal) = $(evaluate.(symXi, tVal))")
get_Ux(L, U, xSym)¶Given a LinearDifferentialOperator L and a VectorBoundaryForm U, constructs the left hand side
$$Ux = M\xi(a) + N\xi(b)$$
of the homogeneous boundary condition $Ux=0$, where
\begin{align*}
Ux &= \begin{bmatrix}
\sum_{j=1}^n (M_{1j}x^{(j-1)}(a) + N_{1j}x^{(j-1)}(b))\\
\vdots\\
\sum_{j=1}^n (M_{mj}x^{(j-1)}(a) + N_{mj}x^{(j-1)}(b))
\end{bmatrix}\\
&= \begin{bmatrix}
M_{11} & \cdots & M_{1n} & N_{11} & \cdots & N_{1n}\\
\vdots & & \vdots & \vdots & & \vdots\\
M_{m1} & \cdots & M_{mn} & N_{m1} & \cdots & N_{mn}
\end{bmatrix} \begin{bmatrix}x(a)\\\vdots\\x^{(n-1)}(a)\\ x(b)\\\vdots\\x^{(n-1)}(b)\end{bmatrix}\\
&= (M:N)\begin{bmatrix}
\xi(a)\\
\xi(b)
\end{bmatrix}.
\end{align*}
function get_Ux(L::LinearDifferentialOperator, U::VectorBoundaryForm, xSym)
(a, b) = L.interval
xi = get_xi(L; symbolic = false, xSym = xSym)
xiEvalA = evaluate.(xi, a)
xiEvalB = evaluate.(xi, b)
M, N = U.M, U.N
Ux = M*xiEvalA + N*xiEvalB
return Ux
end
Parameters
L: LinearDifferentialOperatorU: VectorBoundaryFormsymX*: SymPy.SymReturns
get_boundaryCondition: Array of NumberExample
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
xSym = t^2+2t
Ux = get_Ux(L, U, xSym)
println("U($symX) = $Ux")
check_adjoint(L, U, adjointU, B)¶Given a boundary value problem $$Lx = 0,\quad Ux=0$$ with linear differential operator $L$ and vector boundary form $U$, a candidate adjoint vector boundary form $U^+$, and the matrix $B$ associated with $L$, checks whether the boundary condition $$U^+x = 0$$ is indeed adjoint to the boundary condition $$Ux=0.$$
function check_adjoint(L::LinearDifferentialOperator, U::VectorBoundaryForm, adjointU::VectorBoundaryForm, B::Array)
(a, b) = L.interval
M, N = U.M, U.N
P, Q = (adjointU.M)', (adjointU.N)'
# Avoid InexactError() when taking inv()
BEvalA = convert(Array{Complex}, evaluate.(B, a))
BEvalB = convert(Array{Complex}, evaluate.(B, b))
left = M * inv(BEvalA) * P
right = N * inv(BEvalB) * Q
# println("left = $left")
# println("right = $right")
tol = set_tol(left, right)
return all(i -> isapprox(left[i], right[i]; atol = tol), 1:length(left)) # Can't use == to deterimine equality because left and right are arrays of floats
end
Parameters
L: LinearDifferentialOperatorU: VectorBoundaryFormadjointU: VectorBoundaryFormB: Array of Numberget_B(L).Returns
check_adjoint: Booltrue if adjointU is indeed adjoint to U, orfalse otherwise.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
# M = [3.9 5.4; 1+2*im 2]
# N = [4.7 8.1 + im; 0.5*im 10]
U = VectorBoundaryForm(M, N)
Uc = get_Uc(U)
# Non-symbolic
B = get_B(L)
BHat = get_BHat(L, B)
H = get_H(U, Uc)
J = get_J(BHat, H)
adjointU = get_adjointCand(J)
check_adjoint(L, U, adjointU, B)
get_adjointU(L, U, pDerivMatrix)¶Given a boundary value problem $$Lx = p_0x^{(n)} + p_1x^{(n-1)} + \cdots + p_{n-1}x^{(1)} + p_n x = 0,\quad Ux=0$$ with linear differential operator $L$ and vector boundary form $U$, an $n\times n$ matrix of derivatives $$\begin{bmatrix}p_0 & \cdots & p_0^{(n-1)}\\ \vdots & \ddots & \vdots\\ p_{n-1} & \cdots & p_{n-1}^{(n-1)}\end{bmatrix},$$ construct $U^+$ such that the boundary condition $U^+=0$ is adjoint to the original boundary condition $Ux=0$.
function get_adjointU(L::LinearDifferentialOperator, U::VectorBoundaryForm, pDerivMatrix=get_pDerivMatrix(L))
B = get_B(L; pDerivMatrix = pDerivMatrix)
BHat = get_BHat(L, B)
Uc = get_Uc(U)
H = get_H(U, Uc)
J = get_J(BHat, H)
adjointU = get_adjointCand(J)
if check_adjoint(L, U, adjointU, B)
return adjointU
else
throw(error("Adjoint found not valid"))
end
end
Parameters
L: LinearDifferentialOperatorU: VectorBoundaryFormpDerivMatrix: Array of Function, Number, or SymPy.#get_pDerivMatrix (SymPy.#), orReturns
get_adjointU: VectorBoundaryFormExample
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
# pFunctions = [t->1 t->t+1 t->t^2+t+1]
# L = LinearDifferentialOperator(pFunctions, interval, symL)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
Helps user to find the roots of an exponential polynomial $\Delta(\lambda)$ where $\lambda\in\mathbb{C}$ by visualizing the roots as the intersections of the level curves $\Re(\Delta) = 0$ and $\Im(\Delta) = 0$.
separate_real_imaginary(delta)¶Separates real and imaginary parts of the symbolic expression of $\Delta(\lambda)$, an exponential polynomial in one variable.
separate_real_imaginary_exp(expr)¶Helper function that deals with the case where the toplevel operation is exponentiation.
# although the function body is the same as "power" and "others", this case is isolated because negative exponents, e.g., factor_list(e^(-im*x)), give PolynomialError('a polynomial expected, got exp(-I*x)',), while factor_list(cos(x)) runs normally
function separate_real_imaginary_exp(expr::SymPy.Sym)
result = real(expr) + im*imag(expr)
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) = SymPy.func(sympyExpExpr).Returns
separate_real_imaginary_exp: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = e^(x+im*y)
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_exp(expr)
separate_real_imaginary_power(expr)¶Helper function that deals with the case where the toplevel operation is power.
# we won't be dealing with cases like x^(x^x)
function separate_real_imaginary_power(expr::SymPy.Sym)
result = real(expr) + im*imag(expr)
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) = SymPy.func(sympyPowerExpr).Returns
separate_real_imaginary_power: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = (x+im*y)^2
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_power(expr)
separate_real_imaginary_mult(expr)¶Helper function that deals with the case where the toplevel operation is multiplication.
function separate_real_imaginary_mult(expr::SymPy.Sym)
terms = args(expr)
result = 1
# if the expanded expression contains toplevel multiplication, the individual terms must all be exponentials or powers
for term in terms
# println("term = $term")
# if term is exponential
if SymPy.func(term) == SymPy.func(sympyExpExpr)
termSeparated = separate_real_imaginary_exp(term)
# if term is power (not sure if this case and the case below overlaps)
elseif SymPy.func(term) == SymPy.func(sympyPowerExpr)
termSeparated = separate_real_imaginary_power(term)
# else, further split each product term into indivdual factors (this case also includes the case where term is a number, which would go into the "constant" below)
else
termSeparated = term # term is a number
# (constant, factors) = factor_list(term)
# termSeparated = constant
# # separate each factor into real and imaginary parts and collect the product of separated factors
# for (factor, power) in factors
# factor = factor^power
# termSeparated = termSeparated * (real(factor) + im*imag(factor))
# end
end
# println("termSeparated = $termSeparated")
# collect the product of separated term, i.e., product of separated factors
result = result * termSeparated
end
result = real(result) + im*imag(result)
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) = SymPy.func(sympyMultExpr).Returns
separate_real_imaginary_mult: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = (x+im*y)*2x
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_mult(expr)
separate_real_imaginary_add(expr)¶Helper function that deals with the case where the toplevel operation is addition.
function separate_real_imaginary_add(expr::SymPy.Sym)
x = symbols("x")
# if the expanded expression contains toplevel addition, the individual terms must all be products or symbols
terms = args(expr)
result = 0
# termSeparated = 0 # to avoid undefined error if there is no else (case incomplete)
for term in terms
# println("term = $term")
# if term is a symbol
if SymPy.func(term) == SymPy.func(x)
termSeparated = term
# if term is exponential
elseif SymPy.func(term) == SymPy.func(sympyExpExpr)
termSeparated = separate_real_imaginary_exp(term)
# if term is a power
elseif SymPy.func(term) == SymPy.func(sympyPowerExpr)
termSeparated = separate_real_imaginary_power(term)
# if term is a product
elseif SymPy.func(term) == SymPy.func(sympyMultExpr)
termSeparated = separate_real_imaginary_mult(term)
# if term is a number
else
termSeparated = term
end
# println("termSeparated = $termSeparated")
result = result + termSeparated
end
result = real(result) + im*imag(result)
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) = SymPy.func(sympAddExpr).Returns
separate_real_imaginary_add: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = (x+im*y) + 2x
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_add(expr)
separate_real_imaginary_power_mult_add(expr)¶Helper function that deals with the cases where the toplevel operation is power, addition, or multiplication.
function separate_real_imaginary_power_mult_add(expr::SymPy.Sym)
if SymPy.func(expr) == SymPy.func(sympyPowerExpr)
result = separate_real_imaginary_power(expr)
elseif SymPy.func(expr) == SymPy.func(sympyMultExpr)
result = separate_real_imaginary_mult(expr)
else #if SymPy.func(expr) == SymPy.func(sympyAddExpr)
result = separate_real_imaginary_add(expr)
# else
# result = expr
end
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) = SymPy.func(sympPowerExpr), SymPy.func(sympMultExpr), or SymPy.func(sympAddExpr).Returns
separate_real_imaginary_power_mult_add: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = (x+im*y) + 2x
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_power_mult_add(expr)
separate_real_imaginary_others(expr)¶Helper function that deals with the case where the toplevel operation is not exponentiation, power, multiplication or addition. In this case, expr must be a single term, e.g., x or cos(2x+1), which is a function wrapping around an expression. So we use the other helper functions to separate the expression wrapped in the function and feed it back to the function.
function separate_real_imaginary_others(expr::SymPy.Sym)
# if the expanded expression is neither of the above, it must be a single term, e.g., x or cos(2x+1), which is a function wrapping around an expression; in this case, use the helper function to clean up the expression and feed it back to the function
term = args(expr)[1]
termCleaned = separate_real_imaginary_power_mult_add(term)
result = subs(expr,args(expr)[1],termCleaned)
result = real(result) + im*imag(result)
return result
end
Parameters
expr: SymPy.SymSymPy.func(expr) != SymPy.func(sympPowerExpr), SymPy.func(sympMultExpr), or SymPy.func(sympAddExpr).Returns
separate_real_imaginary_others: SymPy.SymExample
x = symbols("x", real = true)
y = symbols("y", real = true)
expr = cos(x+im*y)
println("func($expr) = $(SymPy.func(expr))")
separate_real_imaginary_others(expr)
separate_real_imaginary(delta)¶The main function that separates the real and imaginary parts of $\Delta$, an exponential polynomial in one variable.
function separate_real_imaginary(delta::SymPy.Sym)
x = symbols("x", real = true)
y = symbols("y", real = true)
freeSymbols = free_symbols(delta)
# check if delta has one and only one free symbol (e.g., global variable lambda)
if length(freeSymbols) == 1
lambda = freeSymbols[1]
# substitute lambda with x+iy
expr = subs(delta, lambda, x+im*y)
# expand the new expression
expr = expand(expr)
if SymPy.func(expr) == SymPy.func(sympyPowerExpr)
# println(expr)
# println("power!")
result = separate_real_imaginary_power(expr)
# println("result = $result")
elseif SymPy.func(expr) == SymPy.func(sympyAddExpr)
# println(expr)
# println("addition!")
result = separate_real_imaginary_add(expr)
# println("result = $result")
elseif SymPy.func(expr) == SymPy.func(sympyMultExpr)
# println(expr)
# println("multiplication!")
result = separate_real_imaginary_mult(expr)
# println("result = $result")
else
# println(expr)
# println("single term!")
result = separate_real_imaginary_others(expr)
# println("result = $result")
end
result = expand(result)
return real(result) + im*imag(result)
else
throw("Delta has more than one variable")
end
end
\begin{align*} \Delta(\lambda)=\lambda + 1 = x + iy + 1. \end{align*}
lambda = symbols("lambda")
delta = lambda + 1
separate_real_imaginary(delta)
\begin{align*} \Delta(\lambda)=e^\lambda = e^{x+iy}=e^x e^{iy} = e^x\cos(y) + ie^x\sin(y)). \end{align*}
delta = e^(lambda)
separate_real_imaginary(delta)
\begin{align*} \Delta(\lambda) &= \lambda^2 = (x+iy)^2 = x^2-y^2 + i2xy. \end{align*}
delta = lambda^2
separate_real_imaginary(delta)
\begin{align*} \Delta(\lambda) &= \cos(\lambda)\\ &= \frac{1}{2}e^{i\lambda} + \frac{1}{2}e^{-i\lambda}\\ &= \frac{1}{2}e^{i(x+iy)} + \frac{1}{2}e^{-i(x+iy)} \\ &= \frac{1}{2}(e^{-y}e^{ix} + e^ye^{-ix})\\ &= \frac{1}{2}e^{-y}(\cos(x) + i\sin(x)) + \frac{1}{2}e^y(\cos(x) - i\sin(x))\\ &= \cos(x)\cosh(y) - i\sin(x)\sinh(y) \end{align*}
delta = cos(lambda)
separate_real_imaginary(delta)
\begin{align*} \Delta(\lambda) &= \cos(\lambda)e^\lambda\\ &= \cos(x+iy)e^{x+iy}\\ &= \cos(x+iy)e^x e^{iy}\\ &= \left[\cos(x)\cosh(y) - i\sin(x)\sinh(y)\right]e^x(\cos(y) + i\sin(y))\\ &= e^x\cos(x)\cosh(y)\cos(y) + e^x\sin(x)\sinh(y)\sin(y) + i(e^x\cos(x)\cosh(y)\sin(y)-e^x\sin(x)\sinh(y)\cos(y)). \end{align*}
delta = cos(lambda)*e^(lambda)
separate_real_imaginary(delta)
\begin{align*} \Delta(\lambda) &= (\lambda^3 + \lambda+2)e^\lambda \end{align*}
delta = (lambda^3+lambda+2)*e^(lambda)
separatedDelta = separate_real_imaginary(delta)
prettyPrint(separatedDelta)
# Verified with WolframAlpha
plot_levelCurves(bivariateDelta; realFunc, imagFunc, xRange, yRange, step, width, height)¶Plots the level curves $\Re(\Delta(x+iy)) = 0$ and $\Im(\Delta(x+iy)) = 0$ in the $\Re-\Im$ plane.
function plot_levelCurves(bivariateDelta::SymPy.Sym; realFunc = real(bivariateDelta), imagFunc = imag(bivariateDelta), xRange = (-INFTY, INFTY), yRange = (-INFTY, INFTY), step = INFTY/1000, width = 1500, height = 1000)
freeSymbols = free_symbols(bivariateDelta)
x = symbols("x", real = true)
y = symbols("y", real = true)
xGridStep = (xRange[2] - xRange[1])/50
yGridStep = (yRange[2] - yRange[1])/50
if freeSymbols == [x, y]
Plots.contour(xRange[1]:step:xRange[2], yRange[1]:step:yRange[2], realFunc, levels=[0], size = (width, height), tickfontsize = 20, seriescolor=:reds, transpose = false, linewidth = 4, linealpha = 1, xticks = xRange[1]:xGridStep:xRange[2], yticks = yRange[1]:yGridStep:yRange[2], grid = true, gridalpha = 0.5)
Plots.contour!(xRange[1]:step:xRange[2], yRange[1]:step:yRange[2], imagFunc, levels=[0], size = (width, height), tickfontsize = 20, seriescolor=:blues, transpose = false, linewidth = 4, linealpha = 1, xticks = xRange[1]:xGridStep:xRange[2], yticks = yRange[1]:yGridStep:yRange[2], grid = true, gridalpha = 0.5)
else
Plots.contour(xRange[1]:step:xRange[2], yRange[1]:step:yRange[2], realFunc, levels=[0], size = (width, height), tickfontsize = 20, seriescolor=:reds, transpose = true, linewidth = 4, linealpha = 1, xticks = xRange[1]:xGridStep:xRange[2], yticks = yRange[1]:yGridStep:yRange[2], grid = true, gridalpha = 0.5)
Plots.contour!(xRange[1]:step:xRange[2], yRange[1]:step:yRange[2], imagFunc, levels=[0], size = (width, height), tickfontsize = 20, seriescolor=:blues, transpose = true, linewidth = 4, linealpha = 1, xticks = xRange[1]:xGridStep:xRange[2], yticks = yRange[1]:yGridStep:yRange[2], grid = true, gridalpha = 0.5)
end
end
Parameters
bivariateDelta: SymPy.SymrealFunc*, imagFunc*: SymPy.Sym or Functionreal(separatedDelta) and imag(separatedDelta). If input manually, they need to be functions in two variables $x$, $y$ where $\lambda = x+iy$.xRange*, yRange*: Tuple{Number, Number}(-INFTY, INFTY).step*: Numberx and y axes, i.e., the step in the sequence of points in x and y to be plotted.width*, height*: NumberReturns
plot_levelCurves: NoneExample
lambda = symbols("lambda")
x = symbols("x", real = true)
y = symbols("y", real = true)
delta = lambda + 1
bivariateDelta = subs(delta, lambda, x+im*y)
# plot_levelCurves(separatedDelta) # somehow causes method error, probably because real(separatedDelta) is a function of x only and imag(separatedDelta) is a function of y only. In this case, we need to input the realFunc and imagFunc manually
plot_levelCurves(bivariateDelta; realFunc = (x, y) -> x + 1, imagFunc = (x, y) -> y)
delta = (lambda^3+lambda+2)*e^(lambda)
bivariateDelta = subs(delta, lambda, x+im*y)
plot_levelCurves(bivariateDelta)
Implement the transform pair (2.15a), (2.15b) on page 10 of "Evolution PDEs and augmented eigenfunctions. Finite interval," given by
\begin{alignat*}{2} F_\lambda: f(x)&\mapsto F(\lambda):\quad F_\lambda(f) &= \begin{cases} F_\lambda^+(f),&\quad\mbox{if $\lambda\in \Gamma_0^+\cup \Gamma_a^+$}\\ F_\lambda^-(f),&\quad\mbox{if $\lambda\in \Gamma_0^-\cup \Gamma_a^-$}\\ \end{cases}\\ f_x: F(\lambda)&\mapsto f(x):\quad f_x(F) &= \int_\Gamma e^{i\lambda x}F(\lambda)\,d\lambda,\quad x\in [0,1]. \end{alignat*}
check_boundaryConditions(L, U, fSym)¶Checks whether $f$ satisfies the boundary conditions, i.e., whether $f\in\Phi$.
function check_boundaryConditions(L::LinearDifferentialOperator, U::VectorBoundaryForm, fSym::Union{SymPy.Sym, Number})
# Checks whether f satisfies the homogeneous boundary conditions
Ux = get_Ux(L, U, fSym)
return check_all(Ux, x->is_approx(x, 0))
end
Parameters
L: LinearDifferentialOperatorU: VectorBoundaryFormfSym: SymPy.Sym or NumberReturns
check_boundaryConditions: Booltrue if $f$ satisfies the homogeneous boundary conditions within a tolerance (TOL) and false otherwise.Example
t = symbols("t")
symPFunctions = [-1 0 0]
interval = (0,1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 0 1]
N = [-1 0; 0 -1]
U = VectorBoundaryForm(M, N)
x = symbols("x")
fSym = sin(2*pi*x)
check_boundaryConditions(L, U, fSym)
get_MPlusMinus(adjointU; symbolic, generic)¶Let $\alpha=e^{2\pi i/n}$; given adjoint vector boundary form associated with two matrices $b^\star$, $\beta^\star$, computes matrices $M^+(\lambda)$, $M^-(\lambda)$ given by \begin{align*} M^+_{kj}(\lambda) &= \sum_{r=0}^{n-1}(-i\alpha^{k-1}\lambda)^rb^\star_{jr}\\ M^-_{kj}(\lambda) &= \sum_{r=0}^{n-1}(-i\alpha^{k-1}\lambda)^r\beta^\star_{jr} \end{align*} as functions of $\lambda$ (for fixed $b^\star$, $\beta^\star$) or their symbolic expressions.
function get_MPlusMinus(adjointU::VectorBoundaryForm; symbolic = false, generic = false)
# these are numeric matrices
bStar, betaStar = adjointU.M, adjointU.N
n = size(bStar)[1]
if symbolic
# return MPlus and MMinus as symbolic expressions with (the global variable) lambda as free variable
lambda = symbols("lambda")
if generic
alpha = symbols("alpha")
else
alpha = e^(2*PI*im/n)
end
MPlusMat = Array{SymPy.Sym}(n,n)
for k = 1:n
for j = 1:n
sumPlus = 0
for r = 0:(n-1)
summandPlus = (-im*alpha^(k-1)*lambda)^r * bStar[j,r+1]
sumPlus += summandPlus
end
sumPlus = simplify(sumPlus)
MPlusMat[k,j] = sumPlus
end
end
MPlusSym = MPlusMat
MMinusMat = Array{SymPy.Sym}(n,n)
for k = 1:n
for j = 1:n
sumMinus = 0
for r = 0:(n-1)
summandMinus = (-im*alpha^(k-1)*lambda)^r * betaStar[j,r+1]
sumMinus += summandMinus
end
sumMinus = simplify(sumMinus)
MMinusMat[k,j] = sumMinus
end
end
MMinusSym = MMinusMat
return (MPlusSym, MMinusSym)
else
alpha = e^(2*pi*im/n)
# if not symbolic, return MPlus and MMinus as functions of lambda
function MPlus(lambda::Number)
MPlusMat = Array{Number}(n,n)
for k = 1:n
for j = 1:n
sumPlus = 0
for r = 0:(n-1)
summandPlus = (-im*alpha^(k-1)*lambda)^r * bStar[j,r+1]
sumPlus += summandPlus
end
MPlusMat[k,j] = sumPlus
end
end
return MPlusMat
end
function MMinus(lambda::Number)
MMinusMat = Array{Number}(n,n)
for k = 1:n
for j = 1:n
sumMinus = 0
for r = 0:(n-1)
summandMinus = (-im*alpha^(k-1)*lambda)^r * betaStar[j,r+1]
sumMinus += summandMinus
end
MMinusMat[k,j] = sumMinus
end
end
return MMinusMat
end
end
return (MPlus, MMinus)
end
Parameters
adjointU: VectorBoundaryFormget_adjointU().symbolic*: Boolgeneric*: Boolsymbolic = true, boolean indicating whether to keep $\alpha=e^{2\pi i/n}$ as a symbol.Returns
get_MPlusMinus: Function or SymPy.Symsymbolic = false and as symbolic expressions if symbolic = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
(MPlusSym, MMinusSym) = get_MPlusMinus(adjointU; symbolic = true, generic = false)
prettyPrint.(MMinusSym)
(MPlus, MMinus) = get_MPlusMinus(adjointU; symbolic = false)
lambda = 1+im
println("MPlus($lambda) = $(MPlus(lambda))")
println("MMinus($lambda) = $(MMinus(lambda))")
println("MPlusSym($lambda) = $(prettyPrint.(evaluate.(MPlusSym, lambda)))")
println("MMinusSym($lambda) = $(prettyPrint.(evaluate.(MMinusSym, lambda)))")
get_M(adjointU; symbolic, generic)¶Computes $M$ given by $$M_{kj}(\lambda) = M^+_{kj}(\lambda) + M^-_{kj}(\lambda)e^{-i\alpha^{k-1}\lambda}$$ as a function of $\lambda$ or its symbolic expression.
function get_M(adjointU::VectorBoundaryForm; symbolic = false, generic = false)
bStar, betaStar = adjointU.M, adjointU.N
n = size(bStar)[1]
if symbolic
# return M as a symbolic expression with lambda as free variable
lambda = symbols("lambda")
if generic
alpha = symbols("alpha")
else
alpha = e^(2*PI*im/n)
end
(MPlusSym, MMinusSym) = get_MPlusMinus(adjointU; symbolic = true, generic = generic)
MLambdaSym = Array{SymPy.Sym}(n,n)
for k = 1:n
for j = 1:n
MLambdaSym[k,j] = simplify(MPlusSym[k,j] + MMinusSym[k,j] * e^(-im*alpha^(k-1)*lambda))
end
end
MSym = simplify.(MLambdaSym)
return MSym
else
alpha = e^(2*pi*im/n)
function M(lambda::Number)
(MPlus, MMinus) = get_MPlusMinus(adjointU)
MPlusLambda, MMinusLambda = MPlus(lambda), MMinus(lambda)
MLambda = Array{Number}(n,n)
for k = 1:n
for j = 1:n
MLambda[k,j] = MPlusLambda[k,j] + MMinusLambda[k,j] * e^(-im*alpha^(k-1)*lambda)
end
end
return MLambda
end
return M
end
end
Parameters
adjointU: VectorBoundaryFormget_adjointU().symbolic*: Boolgeneric*: Boolsymbolic = true, boolean indicating whether to keep $\alpha=e^{2\pi i/n}$ as a symbol.Returns
get_M: Function or SymPy.Symsymbolic = false and as symbolic expression if symbolic = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
MSym = get_M(adjointU; symbolic = true, generic = false)
prettyPrint.(MSym)
M = get_M(adjointU; symbolic = false)
lambda = 1+im
println("M($lambda) = $(M(lambda))")
println("MSym($lambda) = $(evaluate.(sym_to_func.(MSym), lambda))")
get_delta(adjointU; symbolic, generic)¶Computes $\Delta := \det(M)$ as a function of $\lambda$ (for fixed adjointU) or its symbolic expression.
function get_delta(adjointU::VectorBoundaryForm; symbolic = false, generic = false)
if symbolic
MSym = get_M(adjointU; symbolic = true, generic = generic)
deltaSym = simplify(SymPy.det(MSym))
return deltaSym
else
function delta(lambda::Number)
M = get_M(adjointU; symbolic = false)
MLambda = convert(Array{Complex}, M(lambda))
return det(MLambda)
end
return delta
end
end
Parameters
adjointU: VectorBoundaryFormget_adjointU().symbolic*: Boolgeneric*: Boolsymbolic = true, boolean indicating whether to keep $\alpha=e^{2\pi i/n}$ as a symbol.Returns
get_delta: Function or SymPy.Symsymbolic = false and as symbolic expression if symbolic = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
deltaSym = simplify(get_delta(adjointU; symbolic = true, generic = false))
prettyPrint(deltaSym)
delta = get_delta(adjointU; symbolic = false)
lambda = 1+im
println("delta($lambda) = $(delta(lambda))")
println("deltaSym($lambda) = $(evaluate(deltaSym, lambda))")
get_Xlj(adjointU, l, j; symbolic, generic)¶Gets $X_{lj}$, which is the $(n-1)\times (n-1)$ submatrix of $M(\lambda)$ whose $11$-entry is the $(l+1)(j+1)$-entry of $M(\lambda)$, as a function of $\lambda$ (for fixed adjointU, $M$, $l$, $j$), or its symbolic expression.
function get_Xlj(adjointU::VectorBoundaryForm, l::Number, j::Number; symbolic = false, generic = false)
bStar, betaStar = adjointU.M, adjointU.N
n = size(bStar)[1]
if symbolic
MSym = get_M(adjointU; symbolic = true, generic = generic)
MBlockSym = [MSym MSym; MSym MSym]
XljSym = MBlockSym[(l+1):(l+1+n-2), (j+1):(j+1+n-2)]
return XljSym
else
M = get_M(adjointU; symbolic = false)
function Xlj(lambda::Number)
MLambda = M(lambda)
MLambdaBlock = [MLambda MLambda; MLambda MLambda]
XljLambda = MLambdaBlock[(l+1):(l+1+n-2), (j+1):(j+1+n-2)]
return XljLambda
end
return Xlj
end
end
Parameters
adjointU: VectorBoundaryFormget_adjointU().l, j: Intsymbolic*: Boolgeneric*: Boolsymbolic = true, boolean indicating whether to keep $\alpha=e^{2\pi i/n}$ as a symbol.Returns
get_Xlj: Array of Function or SymPy.Symsymbolic = false and as symbolic expressions if symbolic = true.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
l, j = 1, 1
XljSym = prettyPrint.(get_Xlj(adjointU, l, j; symbolic = true))
Xlj = get_Xlj(adjointU, l, j; symbolic = false)
lambda = 1+im
println("Xlj($lambda) = $(Xlj(lambda))")
println("XljSym($lambda) = $(evaluate.(XljSym, lambda))")
get_ChebyshevIntegral(l, f; symbolic, lambda, alpha)¶Computes $$\int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx$$ by approximating $f(x)$ using Chebyshev polynomials and obtaining an explicit expression for the integral.
get_ChebyshevTermIntegral(n; symbolic, c)¶Computes $\tilde{T}_n(c)$ given by \begin{alignat*}{2} \tilde{T}_n(c) = \int_0^\pi e^{-ic\cos\theta}\cos(n\theta)\sin\theta\,d\theta &= \begin{cases} 2 &\quad\mbox{if $n=0$}\\ 0 &\quad\mbox{if $n=1$}\\ \frac{(-1)^{n+1}-1}{n^2-1} &\quad\mbox{if $n\geq 2$} \end{cases}, &\quad\mbox{if $c=0$}\\ &= \sum_{m=1}^{n+1}\alpha(m,n)\left[\frac{e^{i\lambda}}{(i\lambda)^m} + (-1)^{m+n}\frac{e^{-i\lambda}}{(i\lambda)^m}\right]&\quad\mbox{if $c\neq 0$}, \end{alignat*} where \begin{align*} \alpha(m,n) = \begin{cases} (-1)^n&\quad\mbox{if $m=1$}\\ (-1)^{n+1}n^2&\quad\mbox{if $m=2$}\\ (-1)^{n+m-1}2^{m-2}n\sum_{k=1}^{n-m+2}\binom{m+k-3}{k-1}\prod_{j=k}^{m+k-3}(n-j)&\quad\mbox{else} \end{cases} \end{align*} or its symbolic expression.
get_alpha(m, n)¶Computes $\alpha(m,n)$.
function get_alpha(m::Int, n::Int)
if m == 1
result = (-1)^n
elseif m == 2
result = (-1)^(n+1)*n^2
else
sum = 0
for k = 1:(n-m+2)
product = 1
for j = k:(m+k-3)
product *= (n-j)
end
sum += binomial(m+k-3, k-1) * product
end
result = (-1)^(n+m-1)*2^(m-2)*n*sum
end
return result
end
Parameters
m, n: IntReturns
get_alpha: NumberExample
n = 3
println("alpha(1,$n) = $(get_alpha(1, n))")
println("alpha(2,$n) = $(get_alpha(2, n))")
println("alpha(3,$n) = $(get_alpha(3, n))")
function get_ChebyshevTermIntegral(n::Int; symbolic = true, c = symbols("c"))
if symbolic
if c == 0
if n == 0
expr = 2
elseif n == 1
expr = 0
else
expr = ((-1)^(n+1)-1)/(n^2-1)
end
else
expr = 0
for m = 1:(n+1)
summand = get_alpha(m, n) * (e^(im*c)/(im*c)^m + (-1)^(m+n)*e^(-im*c)/(im*c)^m)
expr += summand
end
end
expr = simplify(expr)
return expr
else
function TTilde(c)
if c == 0
if n == 0
result = 2
elseif n == 1
result = 0
else
result = ((-1)^(n+1)-1)/(n^2-1)
end
else
result = 0
for i = 1:(n+1)
summand = get_alpha(i, n) * (e^(im*c)/(im*c)^i + (-1)^(i+n)*e^(-im*c)/(im*c)^i)
result += summand
end
end
return result
end
return TTilde
end
end
Parameters
n: Intsymbolic*: Boolc*: SymPy.Symsymbolic = true, argument of $\tilde{T}_n$. Default to symbols("c").Returns
get_ChebyshevTermintegral: Function or SymPy.Symsymbolic = false and as symbolic expression if symbolic = true.Example
# alpha = e^(2pi*im/1)
# lambda = symbols("lambda")
# c = alpha^(l-1)*lambda/2
chebTermIntSym = get_ChebyshevTermIntegral(2; symbolic = true)
prettyPrint(chebTermIntSym)
chebTermInt = get_ChebyshevTermIntegral(2; symbolic = false)
alpha = e^(2pi*im/1)
lambda = 1+im
l = 2
c = alpha^(l-1)*lambda/2
println("chebTermInt($(prettyPrint(c))) = $(chebTermInt(c))")
println("chebTermIntSym($(prettyPrint(c))) = $(evaluate(chebTermIntSym, c))")
get_ChebyshevCoefficients(f)¶Obtains the coefficients in the Chebyshev approximation of $f$ on $[0,1]$.
Note that the coefficients are obtained by shifting the interval $[0,1]$ to $[-1,1]$ and then back. That is, define $g:[0,1]\to [-1,1]$ by $g(x) = 2x-1$. For $t\in [-1,1]$, define $q(t)=f\circ g^{-1}(t) = f(\frac{t+1}{2}) =: f(x)$ for $x\in [0,1]$. Then the returned coefficients are $\{b_0,\ldots, b_N\}$ in $$f(x) = f\circ g^{-1}(g(x)) = q(g(x)) = q(t) = \sum_{n=0}^N b_n T_n(t) = \sum_{n=0}^N b_n T_n(g(x)) = \sum_{n=0}^N b_n T_n(2x-1)$$ instead of $\{a_0,\ldots, a_N\}$ in $$f(x) = \sum_{n=0}^N a_n T_n(x).$$
function get_ChebyshevCoefficients(f::Union{Function,Number})
fCheb = ApproxFun.Fun(f, 0..1) # Approximate f on [0,1] using chebyshev polynomials
chebCoefficients = ApproxFun.coefficients(fCheb) # get coefficients of the Chebyshev polynomial
return chebCoefficients
end
Parameters
f: Function or NumberReturns
get_ChebyshevCoefficients: Array of NumberExample
f(x) = x^2+1
fChebCoeffs = get_ChebyshevCoefficients(f)
get_ChebyshevIntegral(l, f; symbolic, lambda, alpha)¶Computes $$\int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx$$ by computing \begin{align*} \frac{1}{2e^{ic}}\int_{-1}^1 e^{-ict} q(t)\,dt = \int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx \end{align*} where $c=\frac{\alpha^{l-1}\lambda}{2}$, $q(t) := f\circ g^{-1}(t)$, and \begin{align*} \int_{-1}^1 e^{-ic t}q(t)\,dt &= \int_{-1}^1 e^{-ic t} \sum_{n=0}^N b_n T_n(t)\,dt\\ &= \sum_{n=0}^N b_n \int_{-1}^1 e^{-ic t} T_n(t)\,dt\\ &= \sum_{n=0}^N b_n \int_{-1}^1 e^{-ic t} \cos(n\cdot\cos^{-1}(t))\,dt\\ &= \sum_{n=0}^N b_n \tilde{T}_n(c). \end{align*}
function get_ChebyshevIntegral(l::Number, f::Union{Function, Number}; symbolic = false, lambda = nothing, alpha = nothing)
fChebCoeffs = get_ChebyshevCoefficients(f)
# Replace coefficients too close to 0 by 0
# fChebCoeffs = [if is_approx(x, 0) 0 else x end for x in fChebCoeffs]
if symbolic
lambda = symbols("lambda")
c = alpha^(l-1)*lambda/2
integralSym = 0
for m = 1:length(fChebCoeffs)
fChebCoeff = fChebCoeffs[m]
if is_approx(fChebCoeff, 0)
continue
else
integralSym += fChebCoeffs[m] * get_ChebyshevTermIntegral(m-1; symbolic = true, c = c)
end
end
integralSym = integralSym/(2*e^(im*c))
integralSym = simplify(integralSym)
return integralSym
else
if isa(lambda, Void) || isa(alpha, Void)
throw("lambda, alpha required")
else
c = alpha^(l-1)*lambda/2
integral = 0
for m = 1:length(fChebCoeffs)
fChebCoeff = fChebCoeffs[m]
if is_approx(fChebCoeff, 0)
continue
else
integral += fChebCoeff * get_ChebyshevTermIntegral(m-1; symbolic = false)(c)
end
end
integral = integral/(2*e^(im*c))
return integral
end
end
end
Parameters
l: Intf: Function or Numbersymbolic* : Boollambda*: SymPy.Sym or Numbersymbolic = true and numeric if symbolic = false. Default to symbols("lambda").alpha*: SymPy.Sym or Numbersymbols("alpha").Returns
get_ChebyshevIntegral : Function or SymPy.Symsymbolic = true and as a function of $\lambda$ if symbolic = false.Example
alpha = e^(2pi*im/1)
l = 2
# f(x) = x^2+1
f(x) = sin(2*pi*x)
chebIntSym = simplify(get_ChebyshevIntegral(l, f; symbolic = true, alpha = alpha))
prettyPrint(chebIntSym)
lambda = 1+im
chebInt = get_ChebyshevIntegral(l, f; symbolic = false, lambda = lambda, alpha = alpha)
# Directly compute the integral
g(x) = e^(-im*alpha^(l-1)*lambda*x)
directInt = quadgk(mult_func(g,f), 0, 1)[1]
println("chebIntegral($lambda) = $chebInt")
println("chebIntegralSym($lambda) = $(evaluate.(chebIntSym, lambda))")
println("directIntegral = $directInt")
get_FPlusMinus(adjointU; symbolic, generic)¶Gets $F^+_\lambda$, $F^-_\lambda$ given by \begin{align*} F^+_\lambda(f) &= \frac{1}{2\pi \Delta(\lambda)} \sum_{l=1}^n\sum_{j=1}(-1)^{(n-1)(l+j)}\det X^{lj}(\lambda)M^+_{1j}(\lambda)\int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx,\\ F^-_\lambda(f) &= \frac{-e^{-i\lambda}}{2\pi \Delta(\lambda)} \sum_{l=1}^n\sum_{j=1}(-1)^{(n-1)(l+j)}\det X^{lj}(\lambda)M^-_{1j}(\lambda)\int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx. \end{align*} as functions of $f$ or their symbolic expressions.
function get_FPlusMinus(adjointU::VectorBoundaryForm; symbolic = false, generic = false)
bStar, betaStar = adjointU.M, adjointU.N
n = size(bStar)[1]
if symbolic
lambda = symbols("lambda")
(MPlusSym, MMinusSym) = get_MPlusMinus(adjointU; symbolic = true, generic = generic)
deltaSym = get_delta(adjointU; symbolic = true, generic = generic)
if generic
alpha = symbols("alpha")
c = symbols("c")
FT = SymFunction("FT[f]")(c)
sumPlusSymGeneric = 0
sumMinusSymGeneric = 0
for l = 1:n
summandPlusSymGeneric = 0
summandMinusSymGeneric = 0
for j = 1:n
XljSym = get_Xlj(adjointU, l, j; symbolic = true, generic = true)
integralSymGeneric = subs(FT, c, alpha^(l-1)*lambda)
summandPlusSymGeneric += (-1)^((n-1)*(l+j)) * SymPy.det(XljSym) * MPlusSym[1,j] * integralSymGeneric
summandMinusSymGeneric += (-1)^((n-1)*(l+j)) * SymPy.det(XljSym) * MMinusSym[1,j] * integralSymGeneric
end
sumPlusSymGeneric += summandPlusSymGeneric
sumMinusSymGeneric += summandMinusSymGeneric
end
FPlusSymGeneric = simplify(1/(2*PI*deltaSym)*sumPlusSymGeneric)
FMinusSymGeneric = simplify((-e^(-im*lambda))/(2*PI*deltaSym)*sumMinusSymGeneric)
return (FPlusSymGeneric, FMinusSymGeneric)
else
alpha = e^(2*PI*im/n)
function FPlusSym(f::Union{Function, Number})
sumPlusSym = 0
for l = 1:n
summandPlusSym = 0
for j = 1:n
XljSym = get_Xlj(adjointU, l, j; symbolic = true)
integralSym = get_ChebyshevIntegral(l, f; symbolic = true, lambda = lambda, alpha = alpha)
summandPlusSym += (-1)^((n-1)*(l+j)) * SymPy.det(XljSym) * MPlusSym[1,j] * integralSym
end
sumPlusSym += summandPlusSym
end
return simplify(1/(2*PI*deltaSym)*sumPlusSym)
end
function FMinusSym(f::Union{Function, Number})
sumMinusSym = 0
for l = 1:n
summandMinusSym = 0
for j = 1:n
XljSym = get_Xlj(adjointU, l, j; symbolic = true)
c = alpha^(l-1)*lambda/2
integralSym = get_ChebyshevIntegral(l, f; symbolic = true, lambda = lambda, alpha = alpha)
summandMinusSym += (-1)^((n-1)*(l+j)) * SymPy.det(XljSym) * MMinusSym[1,j] * integralSym
end
sumMinusSym += summandMinusSym
end
return simplify((-e^(-im*lambda))/(2*PI*deltaSym)*sumMinusSym)
end
return (FPlusSym, FMinusSym)
end
else
alpha = e^(2pi*im/n)
(MPlus, MMinus) = get_MPlusMinus(adjointU; symbolic = false)
function FPlus(lambda::Number, f::Union{Function, Number})
MPlusLambda, MMinusLambda = MPlus(lambda), MMinus(lambda)
M = get_M(adjointU; symbolic = false)
MLambda = convert(Array{Complex}, M(lambda))
deltaLambda = det(MLambda) # or deltaLambda = (get_delta(adjointU))(lambda)
sumPlus = 0
for l = 1:n
summandPlus = 0
for j = 1:n
Xlj = get_Xlj(adjointU, l, j; symbolic = false)
XljLambda = convert(Array{Complex}, Xlj(lambda))
integral = get_ChebyshevIntegral(l, f; symbolic = false, lambda = lambda, alpha = alpha)
summandPlus += (-1)^((n-1)*(l+j)) * det(XljLambda) * MPlusLambda[1,j] * integral
end
sumPlus += summandPlus
end
return 1/(2pi*deltaLambda)*sumPlus
end
function FMinus(lambda::Number, f::Union{Function, Number})
MPlusLambda, MMinusLambda = MPlus(lambda), MMinus(lambda)
M = get_M(adjointU; symbolic = false)
MLambda = convert(Array{Complex}, M(lambda))
deltaLambda = det(MLambda) # or deltaLambda = (get_delta(adjointU))(lambda)
sumMinus = 0
for l = 1:n
summandMinus = 0
for j = 1:n
Xlj = get_Xlj(adjointU, l, j)
XljLambda = convert(Array{Complex}, Xlj(lambda))
integral = get_ChebyshevIntegral(l, f; symbolic = false, lambda = lambda, alpha = alpha)
summandMinus += (-1)^((n-1)*(l+j)) * det(XljLambda) * MMinusLambda[1,j] * integral
end
sumMinus += summandMinus
end
return (-e^(-im*lambda))/(2pi*deltaLambda)*sumMinus
end
return (FPlus, FMinus)
end
end
Parameters
adjointU: VectorBoundaryFormsymbolic*: Boolgeneric*: Boolsymbolic = true, boolean indicating whether $f$ and $\alpha=e^{2\pi i/n}$ are kept as generic symbols.Returns
get_FPlusMinusLambda: Tuple of Function or SymPy.Symsymbolic = true (where $\text{FT}[f]$ indicates the Fourier transform integral of $f$, $\int_0^1 e^{-i\alpha^{l-1}\lambda x}f(x)\,dx$) and as functions in $\lambda$ and $f$ if symbolic = false.Example
t = symbols("t")
symPFunctions = [1 t+1 t^2+t+1]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0; 2 0]
N = [0 2; 0 1]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
# generic f
(FPlusSymGeneric, FMinusSymGeneric) = get_FPlusMinus(adjointU; symbolic = true, generic = true)
prettyPrint(FPlusSymGeneric)
f(x) = x^2-1
# Keep lambda as a free symbol
(FPlusSym, FMinusSym) = get_FPlusMinus(adjointU; symbolic = true)
prettyPrint(FPlusSym(f))
lambda = 1+im
(FPlus, FMinus) = get_FPlusMinus(adjointU; symbolic = false)
println("FPlus($lambda, f) = $(FPlus(lambda, f))")
println("FMinus($lambda, f) = $(FMinus(lambda, f))")
println("FPlusSym($lambda, f) = $(evaluate.(FPlusSym(f), lambda))")
println("FMinusSym($lambda, f) = $(evaluate.(FMinusSym(f), lambda))")
t = symbols("t")
symPFunctions = [1 0 0 0]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0 0; 0 0 0; 0 -1/2 0]
N = [0 0 0; 1 0 0; 0 1 0]
U = VectorBoundaryForm(M, N)
adjointU = get_adjointU(L, U)
delta = get_delta(adjointU; symbolic = true, generic = true)
prettyPrint(simplify(delta))
# (FPlusSymGeneric, FMinusSymGeneric) = get_FPlusMinus(adjointU; symbolic = true, generic = true)
# prettyPrint(simplify(FPlusSymGeneric))
get_gamma(a, n, zeroList; infty, nGon)¶Computes the contour $\Gamma$ given by \begin{align*} \Gamma &:= \Gamma_0\cup \Gamma_a, \end{align*} where \begin{align*} \Gamma_a^{\pm} &:= \partial(\{\lambda\in\mathbb{C}^{\pm}:\, \Re(a\lambda^n)>0\}\setminus \bigcup_{\substack{\sigma\in\mathbb{C};\\\Delta(\sigma)=0}} D(\sigma, 2\epsilon))\quad\mbox{($D$ for disk)},\\ \Gamma_a &:= \Gamma_a^+\cup \Gamma_a^-,\\ \Gamma_0^+ &:= \bigcup_{\substack{\sigma\in\overline{\mathbb{C}^+};\\ \Delta(\sigma)=0}} C(\sigma, \epsilon)\quad\mbox{($C$ for circle)},\\ \Gamma_0^- &:= \bigcup_{\substack{\sigma\in\mathbb{C}^-;\\ \Delta(\sigma)=0}} C(\sigma, \epsilon),\\ \Gamma_0 &:= \Gamma_0^+\cup \Gamma_0^-. \end{align*}
get_gammaAAngles(a, n; symbolic)¶Finds the angles (in $[-2\pi, 2\pi)$) representing the lines through the origin that mark the beginning and end of the sectors encircled by $\Gamma_a$, boundary of the domain $\{\lambda\in \mathbb{C}: \Re(a\lambda^n)>0\}$.
Note that we choose the interval $[-2\pi, 2\pi)$ to ensure that "z is in a sector" if and only if "angle(z) is greater than or equal to the start of the sector and smaller than or equal to the end of the sector.
function get_gammaAAngles(a::Number, n::Int; symbolic = false)
# thetaA = argument(a)
thetaA = angle(a)
if symbolic
thetaStartList = Array{SymPy.Sym}(n) # List of angles that characterize where domain sectors start
thetaEndList = Array{SymPy.Sym}(n) # List of angles that characterize where domain sectors end
k = symbols("k")
counter = 0
while (2pi*counter + pi/2 - thetaA)/n < 2pi
# Substituting counter for k
# while SymPy.N(subs((2PI*k + PI/2 - thetaA)/n, k, counter)) < 2pi
thetaStart = (2*PI*counter - PI/2 - rationalize(thetaA/pi)*PI)/n
thetaEnd = (2*PI*counter + PI/2 - rationalize(thetaA/pi)*PI)/n
counter += 1
thetaStartList[counter] = thetaStart
thetaEndList[counter] = thetaEnd
end
else
thetaStartList = Array{Number}(n)
thetaEndList = Array{Number}(n)
k = 0
while (2pi*k + pi/2 - thetaA)/n < 2pi
thetaStart = (2pi*k - pi/2 - thetaA)/n
thetaEnd = (2pi*k + pi/2 - thetaA)/n
k += 1
thetaStartList[k] = thetaStart
thetaEndList[k] = thetaEnd
end
end
return (thetaStartList, thetaEndList)
end
Parameters
a: Numberget_input().n: Intsymbolic*: BoolReturns
get_gammaAAngles: Tuple{Array, Array}Example
n = 2
a = 1
gammaAAnglesSym = get_gammaAAngles(a, n; symbolic = true)
println("n=$n, a=$a, gammaAAnglesSym = $gammaAAnglesSym")
gammaAAngles = get_gammaAAngles(a, n; symbolic = false)
println("n=$n, a=$a, gammaAAngles = $gammaAAngles")
n = 3
a = im
gammaAAnglesSym = get_gammaAAngles(a, n; symbolic = true)
println("n=$n, a=$a, gammaAAnglesSym = $gammaAAnglesSym")
n = 3
a = -im
gammaAAnglesSym = get_gammaAAngles(a, n; symbolic = true)
println("n=$n, a=$a, gammaAAnglesSym = $gammaAAnglesSym")
get_gammaAAnglesSplit(a, n; symbolic)¶On top of the sectors characterized by the array of start angles and the array of end angles obtained via get_gammaAAngles(), splits the sectors that contain the real line.
function get_gammaAAnglesSplit(a::Number, n::Int; symbolic = false)
(thetaStartList, thetaEndList) = get_gammaAAngles(a, n; symbolic = symbolic)
# Split sectors that contain the positive half of the real line (angle = 0)
zeroIndex = find(i -> ((is_approxLess(thetaStartList[i], 0) && is_approxLess(0, thetaEndList[i]))), 1:length(thetaStartList))
if !isempty(zeroIndex)
index = zeroIndex[1]
# Insert 0 after thetaStart
splice!(thetaStartList, (index+1):index, 0)
# Insert 0 before thetaEnd
splice!(thetaEndList, index:(index-1), 0)
end
# Split sectors that contain the negative half of the real line (angle = pi)
piIndex = find(i -> ((is_approxLess(thetaStartList[i], pi) && is_approxLess(pi, thetaEndList[i]))), 1:length(thetaStartList))
if !isempty(piIndex)
index = piIndex[1]
if symbolic
# Insert pi after thetaStart
splice!(thetaStartList, (index+1):index, pi)
# Insert pi before thetaEnd
splice!(thetaEndList, index:(index-1), pi)
else
# Use pi*1 instead of pi to avoid "<= not defined for Irrational{:pi}" error in get_gamma()
splice!(thetaStartList, (index+1):index, pi*1)
splice!(thetaEndList, index:(index-1), pi*1)
end
end
return (thetaStartList, thetaEndList)
end
Parameters
a: Numberget_input().n: Intsymbolic*: BoolReturns
get_gammaAAnglesSplit: Tuple{Array, Array}get_gammaAAngles(), the arrays of angles reflect the fact that sectors containing the real line are splitted into upper half and lower half.Example
n = 2
a = 1
gammaAAnglesSplitSym = get_gammaAAnglesSplit(a, n; symbolic = true)
println("n=$n, a=$a, gammaAAnglesSplitSym = $gammaAAnglesSplitSym")
gammaAAnglesSplit = get_gammaAAnglesSplit(a, n; symbolic = false)
println("n=$n, a=$a, gammaAAnglesSplit = $gammaAAnglesSplit")
n = 3
a = im
gammaAAnglesSplitSym = get_gammaAAnglesSplit(a, n; symbolic = true)
println("n=$n, a=$a, gammaAAnglesSplitSym = $gammaAAnglesSplitSym")
gammaAAnglesSplit = get_gammaAAnglesSplit(a, n; symbolic = false)
println("n=$n, a=$a, gammaAAnglesSplit = $gammaAAnglesSplit")
pointOnSector(z, sectorAngles)¶Determines whether a point $z$ is on the boundary of a sector characterized by a start angle and an end angle.
function pointOnSector(z::Number, sectorAngles::Tuple{Number, Number})
(startAngle, endAngle) = sectorAngles
return is_approx(argument(z), startAngle) || is_approx(argument(z), endAngle) || is_approx(angle(z), startAngle) || is_approx(angle(z), endAngle)
end
Parameters
z: NumbersectorAngles: Tuple{Number, Number}Returns
pointOnSector: Booltrue if z is on the boundary of the sector given by sectorAngles and false otherwise.Example
z = 1+im
sectorAngles = (-pi/4, pi/4)
pointOnSector(z, sectorAngles)
pointInSector(z, sectorAngles)¶Determines whether a point $z$ is in the interior of a sector characterized by a start angle and an end angle.
function pointInSector(z::Number, sectorAngles::Tuple{Number, Number})
(startAngle, endAngle) = sectorAngles
# First check if z is on the sector boundary
if pointOnSector(z, sectorAngles)
return false
else
# angle(z) would work if it's in the sector with positive real parts and both positive and negative imaginary parts; argument(z) would work if it's in the sector with negative real parts and both positive and negative imaginary parts
return (angle(z) > startAngle && angle(z) < endAngle) || (argument(z) > startAngle && argument(z) < endAngle) # no need to use is_approxLess because the case of approximatedly equal is already checked in pointOnSector
end
end
Parameters
z: NumbersectorAngles: Tuple{Number, Number}Returns
pointInSector: Booltrue if z is interior to the sector given by sectorAngles and false otherwise.Example
z = 1+im
sectorAngles = (-pi/4, pi/4)
println("pointInSector($z, $sectorAngles) = $(pointInSector(z, sectorAngles))")
z = 1
println("pointInSector($z, $sectorAngles) = $(pointInSector(z, sectorAngles))")
z = -1
println("pointInSector($z, $sectorAngles) = $(pointInSector(z, sectorAngles))")
pointExSector(z, sectorAngles)¶Determines whether a point $z$ is in the exterior of a sector characterized by a start angle and an end angle.
function pointExSector(z::Number, sectorAngles::Tuple{Number, Number})
return !pointOnSector(z, sectorAngles) && !pointInSector(z, sectorAngles)
end
Parameters
z: NumbersectorAngles: Tuple{Number, Number}Returns
pointExSector: Booltrue if z is exterior to the sector given by sectorAngles and false otherwise.Example
z = 1+im
sectorAngles = (-pi/4, pi/4)
println("pointExSector($z, $sectorAngles) = $(pointExSector(z, sectorAngles))")
z = 1
println("pointExSector($z, $sectorAngles) = $(pointExSector(z, sectorAngles))")
z = -1
println("pointExSector($z, $sectorAngles) = $(pointExSector(z, sectorAngles))")
get_epsilon(zeroList, a, n)¶Given a list of zeroes of $\Delta(\lambda)$, using the arrays of angles that characterize the starts and ends of sectors, computes an appropriate value for the radius $\epsilon$ of the circle to be drawn around each zero of $\Delta(\lambda)$ in the contours $\Gamma_0^{\pm}$. This is the minimum of pairwise distances between zeroes that are not interior to any sector (since interior zeroes would not matter in any way) and the distances from any of these zeroes to any line that mark the boundary of some sector.
function get_epsilon(zeroList::Array, a::Number, n::Int)
(thetaStartList, thetaEndList) = get_gammaAAnglesSplit(a, n; symbolic = false)
thetaStartEndList = collect(Iterators.flatten([thetaStartList, thetaEndList]))
truncZeroList = []
for zero in zeroList
# If zero is interior to any sector, discard it
if any(i -> pointInSector(zero, (thetaStartList[i], thetaEndList[i])), 1:n)
else # If not, append it to truncZeroList
append!(truncZeroList, zero)
end
end
# List of distance between each zero and each line marking the boundary of some sector
pointLineDistances = [get_distancePointLine(z, theta) for z in zeroList for theta in thetaStartEndList]
if length(truncZeroList)>1
# List of distance between every two zeroes
pairwiseDistances = [norm(z1-z2) for z1 in zeroList for z2 in truncZeroList]
else
pairwiseDistances = []
end
distances = collect(Iterators.flatten([pairwiseDistances, pointLineDistances]))
# Distances of nearly 0 could be instances where the zero is actually on some sector boundary
distances = filter(x -> !is_approx(x, 0), distances)
epsilon = minimum(distances)/5
return epsilon
end
Parameters
zeroList: Array of Numberplot_levelCurves().a: Numberget_input().n: IntReturns
get_epsilon: NumberExample
n = 2
a = 1
zeroList = [1+sqrt(3)*im, 2+2*sqrt(3)*im, 0+0*im, 0+5*im, 0-5*im]
get_epsilon(zeroList, a, n)
get_nGonAroundZero(zero, epsilon, n)¶Given a zero of $\Delta(\lambda)$, returns an array of $n$ complex numbers representing the vertices of an $n$-gon around the zero in the complex plane; each vertex is of distance $\epsilon$ from the zero.
function get_nGonAroundZero(zero::Number, epsilon::Number, n::Int)
z = zero
theta = argument(zero)
deltaAngle = 2pi/n
vertices = []
for i = 1:n
newAngle = pi-deltaAngle*(i-1)
vertex = z + epsilon*e^(im*(theta+newAngle))
append!(vertices, vertex)
end
# vertices = vcat(vertices, vertices[1])
return vertices
end
Parameters
zero: Numberepsilon: Numberget_epsilon(). Distance between each vertex of the $n$-gon from zero.n: IntReturns
get_nGonAroundZero: Array of NumberExample
zero = im
epsilon = 1
n = 4
prettyPrint.(get_nGonAroundZero(zero, epsilon, n))
get_gamma(a, n, zeroList; infty, nGon)¶Finds the contours $\Gamma_a^+$, $\Gamma_a^-$, $\Gamma_0^+$, $\Gamma_0^-$ as arrays of points ordered in the order they will be integrated over.
function get_gamma(a::Number, n::Int, zeroList::Array; infty = INFTY, nGon = 8)
(thetaStartList, thetaEndList) = get_gammaAAnglesSplit(a, n; symbolic = false)
nSplit = length(thetaStartList)
gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus = [], [], [], []
epsilon = get_epsilon(zeroList, a, n)
for i in 1:nSplit
thetaStart = thetaStartList[i]
thetaEnd = thetaEndList[i]
# Initialize the boundary of each sector with the ending boundary, the origin, and the starting boundary (start and end boundaries refer to the order in which the boundaries are passed if tracked counterclockwise)
initialPath = [infty*e^(im*thetaEnd), 0+0*im, infty*e^(im*thetaStart)]
initialPath = convert(Array{Complex{Float64}}, initialPath)
if thetaStart >= 0 && thetaStart <= pi && thetaEnd >= 0 && thetaEnd <= pi # if in the upper half plane, push the boundary path to gamma_a+
push!(gammaAPlus, initialPath) # list of lists
else # if in the lower half plane, push the boundary path to gamma_a-
push!(gammaAMinus, initialPath)
end
end
# Sort the zeroList by norm, so that possible zero at the origin comes last. We need to leave the origin in the initial path unchanged until we have finished dealing with all non-origin zeros because we use the origin in the initial path as a reference point to decide where to insert the deformed path
zeroList = sort(zeroList, lt=(x,y)->!isless(norm(x), norm(y)))
for zero in zeroList
# println(zero)
# If zero is not at the origin
if !is_approx(zero, 0+0*im)
# Draw an n-gon around it
vertices = get_nGonAroundZero(zero, epsilon, nGon)
# If zero is on the boundary of some sector
if any(i -> pointOnSector(zero, (thetaStartList[i], thetaEndList[i])), 1:nSplit)
# Find which sector(s) zero is on
indices = find(i -> pointOnSector(zero, (thetaStartList[i], thetaEndList[i])), 1:nSplit)
# If zero is on the boundary of one sector
if length(indices) == 1
# if vertices[2] is interior to any sector, include vertices on the other half of the n-gon in the contour approximation
z0 = vertices[2]
if any(i -> pointInSector(z0, (thetaStartList[i], thetaEndList[i])), 1:nSplit)
# Find which sector vertices[2] is in
index = find(i -> pointInSector(z0, (thetaStartList[i], thetaEndList[i])), 1:nSplit)[1]
else # if vertices[2] is exterior, include vertices on this half of the n-gon in the contour approximation
# Find which sector vertices[length(vertices)] is in
z1 = vertices[length(vertices)]
index = find(i -> pointInSector(z1, (thetaStartList[i], thetaEndList[i])), 1:nSplit)[1]
end
thetaStart = thetaStartList[index]
thetaEnd = thetaEndList[index]
# Find all vertices exterior to or on the boundary of this sector, which would form the nGonPath around the zero
nGonPath = vertices[find(vertex -> !pointInSector(vertex, (thetaStart, thetaEnd)), vertices)]
# If this sector is in the upper half plane, deform gamma_a+
if thetaStart >= 0 && thetaStart <= pi && thetaEnd >= 0 && thetaEnd <= pi
gammaAPlusIndex = find(path -> (is_approx(argument(zero), argument(path[1])) || is_approx(argument(zero), argument(path[length(path)]))), gammaAPlus)[1]
deformedPath = copy(gammaAPlus[gammaAPlusIndex])
if any(i -> is_approx(argument(zero), thetaStartList[i]) || is_approx(angle(zero), thetaStartList[i]), 1:nSplit) # if zero is on the starting boundary, insert the n-gon path after 0+0*im
splice!(deformedPath, length(deformedPath):(length(deformedPath)-1), nGonPath)
else # if zero is on the ending boundary, insert the n-gon path before 0+0*im
splice!(deformedPath, 2:1, nGonPath)
end
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAPlus[gammaAPlusIndex] = deformedPath
else # if sector is in the lower half plane, deform gamma_a-
# # Find all vertices interior to or on the boundary of this sector, which would form the nGonPath around the zero
# nGonPath = vertices[find(vertex -> !pointExSector(vertex, (thetaStart, thetaEnd)), vertices)]
gammaAMinusIndex = find(path -> (is_approx(argument(zero), argument(path[1])) || is_approx(argument(zero), argument(path[length(path)]))), gammaAMinus)[1]
deformedPath = copy(gammaAMinus[gammaAMinusIndex])
if any(i -> is_approx(argument(zero), thetaStartList[i]) || is_approx(angle(zero), thetaStartList[i]), 1:nSplit) # if zero is on the starting boundary, insert the n-gon path after 0+0*im
splice!(deformedPath, length(deformedPath):(length(deformedPath)-1), nGonPath)
else # if zero is on the ending boundary, insert the n-gon path before 0+0*im
splice!(deformedPath, 2:1, nGonPath)
end
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAMinus[gammaAMinusIndex] = deformedPath
end
else # If zero is on the boundary of two sectors, then it must be on the real line, and we need to deform two sectors
# Find out which vertices are in the lower half plane
nGonPath = vertices[find(vertex -> !pointInSector(vertex, (0, pi)), vertices)]
for index in indices
thetaStart = thetaStartList[index]
thetaEnd = thetaEndList[index]
# If this is the sector in the upper half plane, deform gamma_a+
if thetaStart >= 0 && thetaStart <= pi && thetaEnd >= 0 && thetaEnd <= pi
gammaAPlusIndex = find(path -> (is_approx(argument(zero), argument(path[1])) || is_approx(argument(zero), argument(path[length(path)]))), gammaAPlus)[1]
deformedPath = copy(gammaAPlus[gammaAPlusIndex])
if is_approx(argument(zero), argument(deformedPath[length(deformedPath)])) # if zero is on the starting boundary, insert the n-gon path after 0+0*im
splice!(deformedPath, length(deformedPath):(length(deformedPath)-1), nGonPath)
else # if zero is on the ending boundary, insert the n-gon path before 0+0*im
splice!(deformedPath, 2:1, nGonPath)
end
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAPlus[gammaAPlusIndex] = deformedPath
else # If this is the sector in the lower half plane, deform gamma_a-
gammaAMinusIndex = find(path -> (is_approx(argument(zero), argument(path[1])) || is_approx(argument(zero), argument(path[length(path)]))), gammaAMinus)[1]
deformedPath = copy(gammaAMinus[gammaAMinusIndex])
if is_approx(argument(zero), argument(deformedPath[length(deformedPath)])) # if zero is on the starting boundary, insert the n-gon path after 0+0*im
splice!(deformedPath, length(deformedPath):(length(deformedPath)-1), nGonPath)
else # if zero is on the ending boundary, insert the n-gon path before 0+0*im
splice!(deformedPath, 2:1, nGonPath)
end
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAMinus[gammaAMinusIndex] = deformedPath
end
end
end
# Sort each sector's path in the order in which they are integrated over
gammaAs = [gammaAPlus, gammaAMinus]
for j = 1:length(gammaAs)
gammaA = gammaAs[j]
for k = 1:length(gammaA)
inOutPath = gammaA[k]
originIndex = find(x->x==0+0*im, inOutPath)[1]
inwardPath = inOutPath[1:(originIndex-1)]
outwardPath = inOutPath[(originIndex+1):length(inOutPath)]
# Sort the inward path and outward path
if length(inwardPath) > 0
inwardPath = sort(inwardPath, lt=(x,y)->!isless(norm(x), norm(y)))
end
if length(outwardPath) > 0
outwardPath = sort(outwardPath, lt=(x,y)->isless(norm(x), norm(y)))
end
inOutPath = vcat(inwardPath, 0+0*im, outwardPath)
inOutPath = convert(Array{Complex{Float64}}, inOutPath)
gammaA[k] = inOutPath
end
gammaAs[j] = gammaA
end
gammaAPlus, gammaAMinus = gammaAs[1], gammaAs[2]
# If zero is interior to any sector (after splitting by real line), ignore it
# If zero is exterior to the sectors, avoid it
elseif all(i -> pointExSector(zero, (thetaStartList[i], thetaEndList[i])), 1:nSplit)
nGonPath = vcat(vertices, vertices[1]) # counterclockwise
nGonPath = convert(Array{Complex{Float64}}, nGonPath)
# If zero is in the upper half plane, add the n-gon path to gamma_0+
if argument(zero) >= 0 && argument(zero) <= pi
push!(gamma0Plus, nGonPath)
else # If zero is in the lower half plane, add the n-gon path to gamma_0-
push!(gamma0Minus, nGonPath)
end
end
else # If zero is at the origin, we deform all sectors and draw an n-gon around the origin
# deform each sector in gamma_a+
for i = 1:length(gammaAPlus)
deformedPath = gammaAPlus[i]
# find the index of the zero at origin in the sector boundary path
index = find(j -> is_approx(deformedPath[j], 0+0*im), 1:length(deformedPath))
# If the origin is not in the path, then it has already been bypassed
if isempty(index)
else # If not, find its index
index = index[1]
end
# create a path around zero (origin); the origin will not be the first or the last point in any sector boundary because it was initialized to be in the middle, and only insertions are performed. Moreover, the boundary path has already been sorted into the order in which they will be integrated over, so squarePath defined below has deformedPath[index-1], deformedPath[index+1] in the correct order.
squarePath = [epsilon*e^(im*argument(deformedPath[index-1])), epsilon*e^(im*argument(deformedPath[index+1]))]
# replace the zero with the deformed path
deleteat!(deformedPath, index) # delete the origin
splice!(deformedPath, index:(index-1), squarePath) # insert squarePath into where the origin was at
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAPlus[i] = deformedPath
end
# deform each sector in gamma_a-
for i = 1:length(gammaAMinus)
deformedPath = gammaAMinus[i]
index = find(j -> is_approx(deformedPath[j], 0+0*im), 1:length(deformedPath))
if isempty(index)
else
index = index[1]
end
squarePath = [epsilon*e^(im*argument(deformedPath[index-1])), epsilon*e^(im*argument(deformedPath[index+1]))]
deleteat!(deformedPath, index)
splice!(deformedPath, index:(index-1), squarePath)
deformedPath = convert(Array{Complex{Float64}}, deformedPath)
gammaAMinus[i] = deformedPath
end
# Draw an n-gon around the origin and add to gamma_0+
vertices = get_nGonAroundZero(zero, epsilon/2, nGon)
nGonPath = vcat(vertices, vertices[1])
nGonPath = convert(Array{Complex{Float64}}, nGonPath)
push!(gamma0Plus, nGonPath)
end
end
return (gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus)
end
Parameters
a: Numberget_input().n: IntzeroList: Array of Numberplot_levelCurves().infty*: NumberINFTY.nGon*: IntReturns
get_gamma: Tuple{Array, Array, Array, Array}Example
n = 2
a = 1
zeroList = [3+3*sqrt(3)*im, 2+2*sqrt(3)*im, 0+0*im, 0+5*im, 0-5*im]
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
println("gammaAPlus = $gammaAPlus")
println("gammaAMinus = $gammaAMinus")
println("gamma0Plus = $gamma0Plus")
println("gamma0Minus = $gamma0Minus")
plot_contour(gamma; infty)¶Visualize the contours $\Gamma_a^+$, $\Gamma_a^-$, $\Gamma_0^+$, $\Gamma_0^-$.
function plot_contour(gamma::Array; infty = INFTY)
sectorPathList = Array{Any}(length(gamma),1)
for i = 1:length(gamma)
# For each sector path in the gamma contour, plot the points in the path and connect them in the order in which they appear in the path
sectorPath = gamma[i]
# labels = map(string, collect(1:1:length(sectorPath)))
sectorPathList[i] = layer(x = real(sectorPath), y = imag(sectorPath), Geom.line(preserve_order=true))
end
coord = Coord.cartesian(xmin=-infty, xmax=infty, ymin=-infty, ymax=infty, fixed=true)
Gadfly.plot(Guide.xlabel("Re"), Guide.ylabel("Im"), coord, sectorPathList...)
end
Parameters
gamma: Arrayget_gamma().infty*: NumberINFTY.Returns
plot_contour: NoneExample
n = 3
a = im
zeroList = [3+3*sqrt(3)*im, 2+2*sqrt(3)*im, 0+0*im, 0+5*im, 0-5*im, 4]
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
gamma = collect(Iterators.flatten([gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus]))
plot_contour(gamma)
n = 2
a = 1
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
gamma = collect(Iterators.flatten([gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus]))
plot_contour(gamma)
n = 4
a = 1
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
gamma = collect(Iterators.flatten([gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus]))
plot_contour(gamma)
solve_IBVP(L, U, a, zeroList, f; FPlusFunc, FMinusFunc, pDerivMatrix, infty)¶Computes the solution \begin{align*} q(x,t) &= f_x\left(e^{-a\lambda^n t}F_\lambda(f)\right)\\ &= \int_{\Gamma_0^+}e^{i\lambda x}e^{-a\lambda^n t}F_\lambda^+(f)\,d\lambda + \int_{\Gamma_a^+}e^{i\lambda x}e^{-a\lambda^n t}F_\lambda^+(f)\,d\lambda + \int_{\Gamma_0^-}e^{i\lambda x}e^{-a\lambda^n t}F_\lambda^-(f)\,d\lambda + \int_{\Gamma_a^-}e^{i\lambda x}e^{-a\lambda^n t}F_\lambda^-(f)\,d\lambda\\ &= \int_{\Gamma_0^+} + \int_{\Gamma_a^+} e^{ix\lambda - at\lambda^n}F_\lambda^+(f)\,d\lambda + \int_{\Gamma_0^-} + \int_{\Gamma_a^-} e^{ix\lambda - at\lambda^n}F_\lambda^-(f)\,d\lambda. \end{align*}
function solve_IBVP(L::LinearDifferentialOperator, U::VectorBoundaryForm, a::Number, zeroList::Array, f::Function; FPlusFunc = lambda->get_FPlusMinus(adjointU; symbolic = false)[1](lambda, f), FMinusFunc = lambda->get_FPlusMinus(adjointU; symbolic = false)[2](lambda, f), pDerivMatrix = get_pDerivMatrix(L), infty = INFTY)
n = length(L.pFunctions)-1
adjointU = get_adjointU(L, U, pDerivMatrix)
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
function q(x,t)
integrandPlus(lambda) = e^(im*lambda*x)*e^(-a*lambda^n*t) * FPlusFunc(lambda)
integrandMinus(lambda) = e^(im*lambda*x)*e^(-a*lambda^n*t) * FMinusFunc(lambda)
tic()
println("integrandPlus = $(integrandPlus(1+im))")
println("integrandMinus = $(integrandMinus(1+im))")
toc()
# Integrate over individual paths in the Gamma contours
println("gamma0Plus = $gamma0Plus")
tic()
integralGamma0Plus = 0
for path in gamma0Plus
println("path = $path")
if length(path) == 0
path = [im,im]
end
tic()
integralGamma0Plus += quadgk(integrandPlus, path...)[1]
toc()
end
toc()
println("int_0_+ = $integralGamma0Plus")
println("gammaAPlus = $gammaAPlus")
tic()
integralGammaAPlus = 0
for path in gammaAPlus
println("path = $path")
if length(path) == 0
path = [im,im]
end
tic()
integralGammaAPlus += quadgk(integrandPlus, path...)[1]
toc()
end
toc()
println("int_a_+ = $integralGammaAPlus")
println("gamma0Minus = $gamma0Minus")
tic()
integralGamma0Minus = 0
for path in gamma0Minus
println("path = $path")
if length(path) == 0
path = [-im,-im]
end
tic()
integralGamma0Minus += quadgk(integrandMinus, path...)[1]
toc()
end
toc()
println("int_0_- = $integralGamma0Minus")
println("gammaAMinus = $gammaAMinus")
tic()
integralGammaAMinus = 0
for path in gammaAMinus
println("path = $path")
if length(path) == 0
path = [-im,-im]
end
tic()
integralGammaAMinus += quadgk(integrandMinus, path...)[1]
toc()
end
toc()
println("int_a_- = $integralGammaAMinus")
return (integralGamma0Plus + integralGammaAPlus + integralGamma0Minus + integralGammaAMinus)
end
return q
end
Parameters
L: LinearDifferentialOperatorU: VectorBoundaryForma: NumberL in the differential equation.zeroList: Arrayf: FunctionFPlusFunc, FMinusFunc: Functionget_FPlusMinus(). Default to lambda-> get_FPlusMinusLambda2(adjointU; symbolic = false)[1](lambda, f). For better performance, the user may usee the symbolic expressions of $F^+$, $F^-$ to directly define FPlusFunc, FMinusFunc as functions.pDerivMatrix*: ArrayL.pFunctions[i]) implemented as a Function, Number, or SymPy.Sym. Default to the output of get_pDerivMatrix(L).infty*: NumberINFTY.Returns
solve_IBVP: Functionq(x,t) that solves the IBVP.Example
Recall that the IBVP is posed as follows.
Define $n$ linearly independent boundary forms $\{B_j: C\to\mathbb{C}\mid j\in\{1,2,\ldots,n\}\}$ \begin{align*} B_j\phi &= \sum_{k=0}^{n-1} \left(b_{jk}\phi^{(k)}(0) + \beta_{jk}\phi^{(k)}(1)\right),\quad j\in \{1,2,\ldots, n\}. \end{align*} Define \begin{align*} \Phi &= \{\phi\in C: B_j\phi = 0\,\forall j \in \{1,2,\ldots, n\}\}. \end{align*} Let $S:\Phi\to C$ be the linear differential operator $$S\phi(x) = (-i)^n \phi^{(n)}(x).$$
Let $a\in\mathbb{C}$ be a constant.
Consider the IBVP \begin{alignat*}{2} (\partial_t + aS)q(x,t) &= 0,\quad&(x,t)\in (0,1)\times (0,T)\\ q(x,0) = f(x)&\in\Phi,\quad &x\in [0,1]\\ q(\cdot, t) &\in\Phi,\quad &t\in [0,T], \end{alignat*}
Case 1.1
Suppose $n=2$. Then $$S\phi(x)= (-i)^2 \phi^{(2)}(x) = -\phi^{(2)}.$$ Suppose $a=e^{i\theta}$ for $\theta\in [-\frac{\pi}{2}, \frac{\pi}{2}]$.
For $\beta_0, \beta_1\in\hat{\mathbb{C}}$ ($\mathbb{C}$ including $0$ and $\infty$), consider the following boundary conditions $\Phi$: \begin{align*} \varphi(0) + \beta_0\varphi(1) &= 0\\ \varphi'(0) + \beta_1\varphi'(1) &= 0. \end{align*}
We note that in complete form,
Thus, $f(x)\in\Phi$ needs to satisfy \begin{align*} Uf = \begin{bmatrix} 1&0&\beta_0&0\\ 0&1&0&\beta_1 \end{bmatrix} \begin{bmatrix}f(0)\\f'(0)\\f(1)\\f'(1)\end{bmatrix} = 0. \end{align*}
n = 2
beta0 = -1
beta1 = -1
theta = 0
a = e^(im*theta)
t = symbols("t")
symPFunctions = [-1 0 0]
interval = (0,1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
b = [1 0; 0 1]
beta = [beta0 0; 0 beta1]
U = VectorBoundaryForm(b, beta)
f(x) = sin(x*2*pi)
x = symbols("x")
fSym = sin(x*2*PI)
check_boundaryConditions(L, U, fSym)
adjointU = get_adjointU(L, U)
delta = get_delta(adjointU; symbolic = true)
lambda = free_symbols(delta)[1]
x = symbols("x", real = true)
y = symbols("y", real = true)
bivariateDelta = subs(delta, lambda, x+im*y)
plot_levelCurves(bivariateDelta; width = 2500, height = 2000)
zeroList = [0, 2*pi, -2*pi]
(gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus) = get_gamma(a, n, zeroList)
gamma = collect(Iterators.flatten([gammaAPlus, gammaAMinus, gamma0Plus, gamma0Minus]))
plot_contour(gamma)
(FPlusSym, FMinusSym) = get_FPlusMinus(adjointU; symbolic = true)
FPlusSymF = prettyPrint(simplify(FPlusSym(f)))
# FPlusSymF = simplify(FPlusSym(f))
function FPlusSymFunc(lambda)
result = (2*(-e^(im*lambda)+1) * (0.785*lambda^8*e^(im*lambda) - 0.785*lambda^8 - 0.014*im*lambda^7*e^(im*lambda) - 0.014*im*lambda^7 + 31.64*lambda^6*e^(im*lambda) - 31.64*lambda^6 + 19.892*im*lambda^5*e^(im*lambda) + 19.892*im*lambda^5 + 773.28*lambda^4*e^(im*lambda) - 773.28*lambda^4 - 7508.255*im*lambda^3*e^(im*lambda) - 7508.255*im*lambda^3 + 138882.961*lambda^2*e^(im*lambda) - 138882.961*lambda^2 + 743198.707*im*lambda*e^(im*lambda) + 743198.707*im*lambda - 1486397.414*e^(im*lambda) + 1486397.414)*e^(-im*lambda))/(pi*lambda^10*(cos(lambda)-1))
return result
end
FMinusSymF = prettyPrint(FMinusSym(f))
function FMinusSymFunc(lambda)
result = (2*(-e^(im*lambda)+1) * (-0.785*lambda^8*e^(im*lambda) + 0.785*lambda^8 + 0.014*im*lambda^7*e^(im*lambda) + 0.014*im*lambda^7 - 31.64*lambda^6*e^(im*lambda) + 31.64*lambda^6 - 19.892*im*lambda^5*e^(im*lambda) - 19.892*im*lambda^5 - 773.28*lambda^4*e^(im*lambda) + 773.28*lambda^4 + 7508.255*im*lambda^3*e^(im*lambda) + 7508.255*im*lambda^3 - 138882.961*lambda^2*e^(im*lambda) + 138882.961*lambda^2 - 743198.707*im*lambda*e^(im*lambda) - 743198.707*im*lambda + 1486397.414*e^(im*lambda) - 1486397.414) * e^(-2*im*lambda))/(pi*lambda^10*(-cos(lambda)+1))
return result
end
q = solve_IBVP(L, U, a, zeroList, f; FPlusFunc = FPlusSymFunc, FMinusFunc = FMinusSymFunc)
tic()
q(1/2, 1/2)
toc()
# quadgk() is slow around zeroes of Delta (poles)
t = 0.1
# Gadfly.plot(x -> real(q(x,t)), 0, 1)
# Slow if FPlusSymFunc, FMinusSymFunc contain many high powers of lambda (long compilation time of integrandPlus, integrandMinus at each x)
# https://github.com/JuliaLang/julia/issues/19158
Case 1.2
Now, for $b_0, b_1\in\hat{\mathbb{C}}$, consider the following boundary conditions $\Phi$: \begin{align*} \varphi(0) + b_0\varphi'(0) &= 0\\ \varphi(1) + b_1\varphi'(1) &= 0. \end{align*} We note that in complete form,
Case 2.1
Suppose $n=3$. Then $$S\phi(x) = (-i)^3 \phi^{(3)}(x) = i\phi^{(3)}.$$ Suppose $a=\pm i$.
For $\beta_0, \beta_1, \beta_2\in\hat{\mathbb{C}}$ ($\mathbb{C}$ including $0$ and $\infty$), consider the following boundary conditions $\Phi$: \begin{align*} \phi(0) + \beta_0\phi(1) &= 0\\ \phi^{(1)}(0) + \beta_1\phi^{(1)}(1) &= 0\\ \phi^{(2)}(0) + \beta_2\phi^{(2)}(1) &= 0. \end{align*}
We note that in complete form,
Case 2.2
Now, for $b_0\in\hat{\mathbb{C}}$, consider the following boundary conditions $\Phi$: \begin{align*} \varphi(0) &= 0\\ \varphi(1) &= 1\\ \varphi^{(1)}(0) + b_0\varphi^{(1)}(1) &= 0. \end{align*} We note that in complete form,
\begin{alignat*}{3} q_t(x,t) + q_{xxx}(x,t) &= 0,\quad &(x,t)&\in (0,1)\times (0,T)\\ q(x,0) &= f(x),\quad &x&\in [0,1]\\ q(0,t) &= 0, \quad &t&\in [0,T]\\ q(1,t) &= 0\quad &t&\in [0,T]\\ q_x(1,t) &= \frac{1}{2}q_x(0,t)\quad &t&\in [0,T]. \end{alignat*}
\begin{align*} F_\lambda^+(f) &= \frac{1}{2\pi\Delta(\lambda)}\left[\text{FT}[f](\lambda)(e^{i\lambda} + 2\alpha e^{-i\alpha\lambda} + 2\alpha^2 e^{-i\alpha^2\lambda}) + \text{FT}[f](\alpha\lambda)(\alpha e^{i\alpha\lambda} - 2\alpha e^{-i\lambda}) \right.\\ &\qquad \left. + \text{FT}[f](\alpha^2\lambda)(\alpha^2e^{i\alpha\lambda} - 2\alpha^2e^{-i\lambda})\right]\\ F_\lambda^-(f) &= \frac{e^{-i\lambda}}{2\pi\Delta(\lambda)}\left[-\text{FT}[f](\lambda)(2+\alpha^2 e^{-\alpha\lambda} + \alpha e^{-\alpha^2\lambda}) - \alpha\text{FT}[f](\alpha\lambda)(2-e^{-i\alpha^2\lambda}) \right.\\ &\qquad \left. - \alpha^2\text{FT}[f](\alpha^2\lambda)(2-e^{-i\alpha\lambda})\right], \end{align*} where \begin{align*} \Delta(\lambda) = e^{i\lambda} + \alpha e^{i\alpha\lambda} + \alpha^2 e^{i\alpha^2\lambda} + 2(e^{-i\lambda} + \alpha e^{-i\alpha\lambda} + \alpha^2e^{-i\alpha^2\lambda}). \end{align*}
n = 3
t = symbols("t")
symPFunctions = [1 0 0 0]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0 0; 0 0 0; 0 1 0]
N = [0 0 0; 1 0 0; 0 -2 0]
U = VectorBoundaryForm(M, N)
c = symbols("c")
FT = SymFunction("FT[f]")(c)
FTc = symbols("FT[f](c)")
alpha = symbols("alpha")
lambda = symbols("lambda")
deltaFormula = e^(im*lambda) + alpha*e^(im*alpha*lambda) + (alpha)^2*e^(im*(alpha)^2*lambda) + 2*(e^(-im*lambda) + alpha*e^(-im*alpha*lambda) + (alpha)^2*e^(-im*(alpha)^2*lambda))
FPlusFormula = 1/(2*PI*deltaFormula) * (FT(lambda)*(e^(im*lambda) + 2*alpha*e^(-im*alpha*lambda) + 2*alpha^2*e^(-im*alpha^2*lambda)) + FT(alpha*lambda)*(alpha*e^(im*alpha*lambda) - 2*alpha*e^(-im*lambda)) + FT(alpha^2*lambda)*(alpha^2*e^(im*alpha^2*lambda) - 2*alpha^2*e^(-im*lambda)))
FMinusFormula = (e^(-im*lambda))/(2*PI*deltaFormula) * (-FT(lambda)*(2+alpha^2*e^(-im*alpha*lambda)+alpha*e^(-im*alpha^2*lambda)) - alpha*FT(alpha*lambda)*(2-e^(-im*alpha^2*lambda)) - alpha^2*FT(alpha^2*lambda)*(2-e^(-im*alpha*lambda)))
adjointU = get_adjointU(L, U)
(FPlusSymGeneric, FMinusSymGeneric) = get_FPlusMinus(adjointU; symbolic = true, generic = true)
prettyPrint(simplify(FPlusSymGeneric))
prettyPrint(simplify(FMinusSymGeneric))
function test_formula(formula1, formula2)
results = [true]
for i = 1:100
println("i = $i")
realLambda = rand(Uniform(-50.0,50.0), 1, 1)[1]
imagLambda = rand(Uniform(-50.0,50.0), 1, 1)[1]
lambdaVal = realLambda + im*imagLambda
println("lambda = $lambdaVal")
println("-----------------------------------------------------------------")
println(FT(lambda))
output1 = SymPy.N(subs(formula1, (FT(lambda), 1), (FT(alpha*lambda), 0), (FT(alpha^2*lambda), 0), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
output2 = SymPy.N(subs(formula2, (FT(lambda), 1), (FT(alpha*lambda), 0), (FT(alpha^2*lambda), 0), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
println("Formula1 = $output1")
println("Formula2 = $output2")
println("-----------------------------------------------------------------")
result = is_approx(output1-output2, 0)
append!(results, result)
println(FT(alpha*lambda))
output1 = SymPy.N(subs(formula1, (FT(lambda), 0), (FT(alpha*lambda), 1), (FT(alpha^2*lambda), 0), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
output2 = SymPy.N(subs(formula2, (FT(lambda), 0), (FT(alpha*lambda), 1), (FT(alpha^2*lambda), 0), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
println("Formula1 = $output1")
println("Formula2 = $output2")
println("-----------------------------------------------------------------")
result = is_approx(output1-output2, 0)
append!(results, result)
println(FT(alpha^2*lambda))
output1 = SymPy.N(subs(formula1, (FT(lambda), 0), (FT(alpha*lambda), 0), (FT(alpha^2*lambda), 1), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
output2 = SymPy.N(subs(formula2, (FT(lambda), 0), (FT(alpha*lambda), 0), (FT(alpha^2*lambda), 1), (alpha, e^(2*pi*im/n)), (lambda, lambdaVal)))
println("Formula1 = $output1")
println("Formula2 = $output2")
println("-----------------------------------------------------------------")
result = is_approx(output1-output2, 0)
append!(results, result)
println("=================================================================")
end
return all(results)
end
test_formula(FPlusFormula, FPlusSymGeneric)
test_formula(FMinusFormula, FMinusSymGeneric)
\begin{alignat*}{3} q_t(x,t) + q_{xxx}(x,t) &= 0,\quad &(x,t)&\in (0,1)\times (0,T)\\ q(x,0) &= f(x),\quad &x&\in [0,1]\\ q(0,t) &= 0, \quad &t&\in [0,T]\\ q(1,t) &= 0\quad &t&\in [0,T]\\ q_x(1,t) &= 0 \quad &t&\in [0,T]. \end{alignat*}
\begin{align*} F_\lambda^+(f) &= \frac{1}{2\pi\Delta(\lambda)}\left[\text{FT}[f](\lambda)(\alpha e^{-\alpha\lambda} + \alpha^2 e^{-i\alpha^2\lambda}) - (\alpha\text{FT}[f](\alpha\lambda) + \alpha^2\text{FT}[f](\alpha^2\lambda))e^{-i\lambda}\right]\\ F_\lambda^-(f) &= \frac{e^{-i\lambda}}{2\pi\Delta(\lambda)}\left[-\text{FT}[f](\lambda) - \alpha\text{FT}[f](\alpha\lambda) - \alpha^2\text{FT}[f](\alpha^2\lambda)\right], \end{align*} where \begin{align*} \Delta(\lambda) = e^{-i\lambda} + \alpha e^{-i\alpha\lambda} + \alpha^2e^{-i\alpha^2\lambda}. \end{align*}
n = 3
t = symbols("t")
symPFunctions = [1 0 0 0]
interval = (0, 1)
symL = SymLinearDifferentialOperator(symPFunctions, interval, t)
L = get_L(symL)
M = [1 0 0; 0 0 0; 0 0 0]
N = [0 0 0; 1 0 0; 0 1 0]
U = VectorBoundaryForm(M, N)
c = symbols("c")
FT = SymFunction("FT[f]")(c)
FTc = symbols("FT[f](c)")
alpha = symbols("alpha")
lambda = symbols("lambda")
deltaFormula = e^(-im*lambda) + alpha*e^(-im*alpha*lambda) + alpha^2*e^(-im*alpha^2*lambda)
FPlusFormula = 1/(2*PI*deltaFormula) * (FT(lambda)*(alpha*e^(-im*alpha*lambda) + alpha^2*e^(-im*alpha^2*lambda)) - (alpha*FT(alpha*lambda) + alpha^2*FT(alpha^2*lambda))*e^(-im*lambda))
FMinusFormula = (e^(-im*lambda))/(2*PI*deltaFormula) * (-FT(lambda) - alpha*FT(alpha*lambda) - alpha^2*FT(alpha^2*lambda))
adjointU = get_adjointU(L, U)
(FPlusSymGeneric, FMinusSymGeneric) = get_FPlusMinus(adjointU; symbolic = true, generic = true)
prettyPrint(simplify(FPlusSymGeneric))
prettyPrint(simplify(FMinusSymGeneric))
test_formula(FPlusFormula, FPlusSymGeneric)
test_formula(FMinusFormula, FMinusSymGeneric)